rfproouced  at  government  expense 


fcFOSR-TR-  G  *  •  1  102 


RESPONSE  OF  SATURATED 
POROUS  NONLINEAR  MATERIALS 
TO  DYNAMIC  LOADINGS 


31  May  1984 


Prepared  for 

Air  Force  Office  of  Scientific  Research 
Washington,  D.C. 


Prepared  under 

Contract  No.  F49620-81-C-0014 


DTIC 

Kwang  Jin  Kirn 
Scott  E.  Blouin 

D 


Applied  Research  Associates,  Inc. 
New  England  Division 
South  Royal ton,  Vermont  05C C° 


-UNCLASSIFIED. 


SECURITY  CLASSIFICATION  OF  this  page 


REPORT  DOCUMENTATION  PAGE 


1*  REPORT  SECURITY  CLASSIFICATION 

UNCLASSIFIED 


ID.  RESTRICTIVE  MARKINGS 


2*.  SECURITY  CLASSIFICATION  AUTHORITY 


2b  OECLASSIFICATION/OOWNGRAOING  schedule 


3  OISTRIOUTION/AVAILABILITY  of  report 

Approved  for  Public  Release; 
Distribution  Unlimited. 


*  performing  organization  report  NUMBERIS) 


5.  MONITORING  ORGANIZATION  REPORT  NUMBERIS) 

AFOSR-TR.  8  4-1  .102 


6*.  NAME  OF  PERFORMING  ORGANIZATION 

applied  research  associates,  in 


)Sb  Office  SYMBOL 
i If  applicable  i 


7a.  NAME  OF  MONITORING  ORGANIZATION 


,  A 


Hr-  OotL 


Sc.  ADDRESS  (City.  Stott  and  ZIP  Codtl 

NEW  ENGLAND  DIVISION 
SOUTH  R0YALT0N,  VT  05068 


7b.  AOORESS  (City.  Slate  and  Zp>  Code) 

/O 


"NS 

).  PPtOtjj 


0  ~Vh?v'---  V 


V.) 


-r- — 

riON  NU 


l 

l 


84  NAME  OF  FUNOING/SPONSORING 

organization  Air  FORCE 
OFFICE  OF  SCIENTIFIC  RESEARCH 


8b  OFFICE  SYMBOL 
(It  applicable) 


AFOSR/NA 


9.  PRot/llREMENT  INSTR 

F49620-81-C-0014 


¥  ID6NTIFICAT 


MBER 


8c  AOORESS  (City,  State  and  7-IP  Code} 

BOLLING  AFB,  DC  20332 


10.  SOURCE  OF  FUNOINC  NOS. 


PROGRAM 
element  NO 

4ll0^p 


11  TITLE  ilnelude  Security  Clauificationi  .  .  . 

RESPONSE  OF  SATURATED  POROUS  NONLINEAR  MATERlAlifiN CL  ASSl  F I 

nVMAUTP  in AT\TKir*c  I 


PROJECT 

NO. 

ED) 


TASK 

NO 


WORK  UNIT 
NO. 


12  personal  autmoais) 

KWANG  J.  KIM,  SCOTT  E. 

BLOUIN 

— — ^  - 

13a  type  of  report 

13b  TIME  COVERED 

14  DATE  OF  REPORT  f Yr  ,  Mo.,  Day) 

15  PAGE  COUNT 

FINAL 

FROM  TO 

31  MAY  1984 

100 

16  SUPPLEMENTARY  NOTATION 

1 17  /  COSATl  COOES  | 

HELD 

GROUP 

SUB.  GR  ^  - 

18.  SUBJECT  TERMS  (Continue  on  revert*  if  necessary  and  identify  by  block  number ; 

'•  Finite  element  modeling,  soil  mechanics,  liquefaction* 
blast-induced  liquefaction, Asoil  dynamics., 


19.  ABSTRACT  (Continue  on  reverse  if  necettary  and  identify  by  block  number } 

y  Past  theoretical  treatments  of  two  phase  media  are  reviewed  and  incorporated  into  the 
two  phase  dynamic  finite  element  code  TPDAP .  A  study  cf  the  response  of  two  phase  porous 
elastic  media  to  dynamic  uniaxial  loadings  is  conducted  and  the  influence  of  loading  shape, 
material  properties,  and  numerical  techniques  and  parameters  are  evaluated.  A  series  of 
calculations  simulating  uniaxial  loadings  of  saturated  soils  having  hysf.eretic  skeletons 
is  conducted.  Zones  of  liquefaction  occuring  in  these  calculations  are  identified  and 
the  mechanics  of  liquefaction  and  influence  of  parametric  variations  on  liquefaction  are 
investigated  and  described.  '  . 


20  OISTRIBUTION/AVAILABILITY  of  abstract 

UN  CL  ASS  I  FIEO/UNL'MlT  C  7?  SAME  “S  APT  G  C/Tlc  USERS  □ 


21  ABSTRACT  SECURITY  CLASSIFICATION 

UNCLASSIFIED 


22*.  NAME  OF  RESPONSIBLE  INDIVIDUAL 

Li.  Col  Lawrence  I)  llnhari.snn 


22b  TELEPHONE  NUMBER 
(Include  .Irco  Code  i 


20. 


*  o  7  —  "i  2  ."s  :'•» 


22c  OFFICE  SYMBOL 


CD  FORM  1473(83  APR 


COITION  OF  1  JAN  73  IS  OBSOLETE 


SE  CUR  IT  V  CLASS  IF  ICAT  1QN  OF  THIS  PAGE 


TABLE  OF  CONTENTS 


SECTION 

PAGE 

1 

INTRODUCTION 

2 

2 

THEORETICAL  DEVELOPMENT  AND  IM¬ 
PLEMENTATION 

5 

3 

TWO  PHASE  DYNAMIC  ANALYSIS  - 
SATURATED  POROUS  ELASTIC  SKELETON 

33 

4 

TWO  PHASE  DYNAMIC  ANALYSIS  - 
SATURATED  POROUS  INELASTIC  SKELETON 

70 

5 

SUMMARY  AND  CONCLUSIONS 

95 

REFERENCES 

99 

SECTION  1 


INTRODUCTION 


Blast  induced  liquefaction  has  been  identified  as  an 
important  mechanism  in  the  development  of  craters  and  ground 
motions  from  large  explosions  in  or  on  saturated  soil  deposits. 
Liquefaction  may  have  dominated  crater  formation  processes  in  the 
saturated  geologies  of  the  Pacific  Proving  Grounds  where  all  the 
high  yield  U.S.  nuclear  cratering  data  were  obtained.  Based  on 
analysis  of  these  and  other  high  explosive  test  data  from  si'-.es 
having  shallow  water  tables  within  the  continental  United  States 
and  Canada,  Blouin  and  Shinn  (1983)  hypothesized  liquefaction  and 
cratering  mechanisms  controlled  by  the  liquefaction  which 
explain  the  unusually  flat  broad  shapes  of  these  craters.  Based 
on  their  hypothesis,  it  was  concluded  that  current  cratering 
prediction  procedures  founded  on  the  Pacific  data  are  not  valid 
for  dry  continental  sites.  The  data  analysis  indicates  that 
water  table  will  influence  the  crater  formation  processes 
whenever  the  scaled  depth  of  the  water  table  is  less  than  about  7 
f t/ton^^ ,  where  the  yield  is  expressed  as  an  equivalent  high 
explosive  half  buried  yield.  It  was  also  concluded,  however, 
that  many  so  called  dry  sites  may  actually  behave  like  wet  sites, 
because  scaled  water  table  depths  for  the  high  yield  nuclear 
threats  of  interest  are  quite  shallow. 

This  previous  work  also  included  formulation  of 
liquefaction  mechanisms  due  to  explosive  loadings.  Application  of 
those  mechanisms  was  limited  by  the  simplifying  assumptions  u«ed 
to  derive  the  analytic  equations. 
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More  precise  methods  to  analyze  and  predict  the 
occurrence  of  blast  induced  liquefaction  must  rely  on  two  phase 
dynamic  numerical  computer  codes.  Such  codes  can  incorporate 
complicated  models  of  real  world  material  properties  with 
accurate  characterization  of  actual  explosive  loacing  functions 
to  determine  the  dynamic  response  of  saturated  geologic  media. 

The  complexity  of  the  real  world  material  properties,  loading 
functions  and  problem  geometries  puts  realistic  depiction  of  the 
dynamic  response  beyond  the  reach  of  conventional  analytic 
techniques . 

In  general  terms,  two  phase  dynamic  codes  compute  the 
response  of  saturated  porous  media  to  static  or  dynamic  loadings. 
The  media  are  modeled  by  combining  the  response  characteristics 
of  the  soil  skeleton  with  the  response  characteristics  of  the 
pore  water  using  two  equations.  The  first  of  these  describes  the 
motion  of  the  bulk  soil-water  mixture  as  a  function  of  the  applied 
load  and  the  second  describes  the  motion  of  the  pore  fluid 
relative  to  the  soil  skeleton  under  the  applied  load.  The 
dynamic  response  of  the  saturated  material  is  expressed  in  terms 
of  the  intergranular  or  effective  stress  in  the  soil  skeleton,  the 
absolute  motion  of  the  soil  skeleton,  the  pressure  in  the  pore 
water,  and  the  motion  of  the  pore  water  relative  to  the  soil 
skeleton.  Presentation  of  the  calculational  results  in  terms  of 
the  response  of  each  separate  phase  makes  the  two  phase  codes 
extremely  powerful  tools  for  the  study  of  blast  induced 
1 iquef act  ion. 

The  objectives  of  this  study  are  to  review  and  condense 
past  theoretical  treatments  of  two  phase  media,  to  incorporate 
these  results  in  the  Two  Phase  Dynamic  Analysis  Program  ( TPDAP ) , 
to  use  TPDAP  to  study  the  response  of  saturated  porous  media  to 
dynamic  loadings,  and  to  identify  and  analyze  liquefaction 
occurring  in  these  calculations.  Procedures  used  in,  and  results 
of  this  study  are  described  in  the  following  three  sections.  The 
objectives  and  content  of  each  section  are  summarized  below. 
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SECTION  2:  Fundamental  theoretical  treatments  of  two  phase 
porous  media  by  Biot  and  others  are  reviewed.  Biot's  fluid 
equation  describing  the  interactions  between  the  pore  fluid  and 
soil  skeleton  is  generalized  and  incorporated  into  the  finite 
element  code  TP DAP .  Also,  a  general  description  of  TPDAP  is 
presented. 

SECTION  3:  A  series  of  parametric  calculations  of  the  response 
of  two  phase  porous  elastic  media  to  dynamic  uniaxial  loadings 
is  conducted.  The  influence  of  various  loadings,  material 
properties  and  numerical  techniques  are  described  and  analyzed. 

SECTION  4:  A  series  of  parametric  calculations  on  saturated 
materials  having  hysteretic  soil  skeletons  is  conducted.  These 
include  a  bilinear  material  skeleton  and  an  actual  sand  skeleton 
Zones  of  liquefaction  occurring  in  these  calculations  are 
identified  and  the  mechanics  of  liquefaction  and  influence  of 
parametric  variations  on  liquefaction  are  described. 


SECTION  2 


THEORETICAL  DEVELOPMENT  AND  IMPLEMENTATION 


BACKGROUND 


The  fundamental  analytic  work  describing  the  behavior  of 
saturated  porous  media  was  performed  by  Biot.  Biot's  results 
were  reported  in  a  series  of  papers  extending  over  many  years 
(e.g.  1956,  1962a  and  1962b).  Other  investigators  have  applied 
Biot's  analytic  results  using  techniques  which  approximate  his 
equations  with  varying  degrees  of  accuracy  and  sophistication 
(e.g.  Ghaboussi  and  Wilson  1972,  Mengi  and  McNiven,  1977).  In 
this  section,  we  summarize  the  development  of  the  two  phase 
equations  of  motion  based  on  Biot's  work  and  their  implementation 
in  the  two  phase  finite  element  code  TPDAP. 


ANALYTIC  DEVELOPMENT 

The  behavior  of  saturated  porous  material  can  be 
described  in  terms  of  two  equations  of  motion;  the  first  is  the 
equation  of  motion  of  the  bulk  soil-water  mixture  and  the  second 
is  the  equation  of  motion  of  the  pore  fluid  within  the  soil 
ske leton . 


Governing  Equation  for  Bulk  Mixture 


The  differential  equation  of  motion  governing  the  bulk 
mixture  is  given  by 
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(2-1) 


o 


13*3 


(1  -  n)ps  tt  +  n  pf  a 


o. •  •  is  the  total  stress  gradient  applied  to  an  infinitesimal 
1-3*3 

element  of  saturated  material  at  some  given  time.  a.  .  .  is 

expressed  in  tensor  notation  and  represents  the  stress  gradient 
in  each  of  3  mutually  perpendicular  coordinates  (e.g.  see 
Mendelson,  1968).  For  instance,  in  the  x  direction, 


o  .  . 

XJ  *3 


+ 


(1  -  n)  pgux  +  n  pf  i3x  (2-2) 


The  term  (1  -  n)pg  is  the  mass  of  the  soil  skeleton  per  unit 
volume  of  saturated  material,  where  n  is  the  porosity  and  Pg  is 
the  mass  density  of  the  solid  grains.  is  the  displacement  of 

the  skeleton  in  the  i  direction  and  vk  is  the  acceleration  of  the 
skeleton  in  the  i  direction.  The  term  npf  is  the  mass  of  pore 
fluid  per  unit  volume  of  saturated  material  where  Pf  is  the  mass 
density  of  the  pore  fluid.  IK  is  the  absolute  displacement  of  the 
pore  fluid  in  the  i  direction  and  IK  the  absolute  acceleration  of 
the  pore  fluid  in  the  i  direction. 

The  bulk  mass  density  of  the  saturated  material,  p,  is 

given  by 


p  =  (1  -  n) ps  +  npf  ( 2-3) 

Substitution  of  the  value  for  (1  -  n)pg  from  Equation  2-3  into 
Equation  2-1  gives 


°ij*j  “  (  P  “  npf^i  +  npf^i  (2_4) 

A  term  is  introduced  which  is  the  apparent  fluid  displacement 
in  the  i  direction  relative  to  the  soil  skeleton  and  is  given  by 
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(2-5) 


w.  =  n (U.  -  u. ) 

1  i  i 

In  seepage  problems,  w^  is  referred  to  as  the  discharge 
displacement.  It  describes  the  discharge  of  fluid  through  a  soil 
mass  of  unit  area.  The  discharge  velocity,  or  apparent  relative 
velocity,  w^,  between  the  soil  particles  and  pore  water  is  the 
velocity  of  water  in  a  discharge  duct  of  unit  area  needed  to 
maintain  the  actual  relative  velocity  in  the  porous  soil  of  the 
same  unit  area.  The  actual  relative  velocity  between  the 
skeleton  and  the  pore  water  is  given  by  w^/n.  Finally,  w  ^  is 
the  apparent  relative  acceleration  between  the  soil  skeleton  and 
pore  water  given  by 


w±  =  nflh  -  ii. )  (  2-6) 

Equation  2-4  can  be  expressed  in  terms  of  the  apparent  relative 
fluid  acceleration  as  simply 


IDrD 


=  pUi 


+  p,w 


f  i 


(2-7) 


Governing  Equation  for  Pore  Fluid 


Biot's  Equations 

The  equation  of  motion  governing  the  pore  fluid  is 
derived  from  Biot's  work  (see  previous  references).  Biot 
expressed  the  pore  pressure  gradient,  tt,^,  as 

n,.  —  p  _  U.  +  D.  (  2-8 ) 

l  l  l 

In  this  equation,  the  pore  pressure  gradient  is  expressed  in 
tensor  notation.  For  example,  the  gradient  in  the  x  direction 
is  given  by 
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(2-9) 


U  +  D 
X  X 


The  term  PfUi  is  the  inertial  force  per  unit  volume  of  pore 
fluid.  represents  the  viscous  friction  force  between  the  pore 

fluid  and  the  soil  skeleton  per  unit  volume  of  pore  fluid. 

Solving  liquation  2-6  for  and  substitution  into  Equation  2-8 
gives 


71 

'i 


—  w .  +  p  ,  U  .  4  D . 

n  l  l  l 


(2-10) 


Biot  showed  that  the  viscous  friction  term,  D  ,  is  a 
function  of  the  excitation  frequency,  w,  the  pore  geometry,  the 
dynamic  viscosity,  u ,  and  the  apparent  relative  velocity  between 
the  pore  fluid  and  the  skeleton,  w-.  in  an  actual  soil  the  flow 
of  pore  water  would  follow  very  complicated  paths  which  are 
difficult  to  describe.  These  flow  paths  would  involve  numerous 
variations  in  direction  and  in  cross  secticnal  area.  Riot 
employed  models  of  the  flow  paths  which  are  gross  Simplifications 
of  the  actual  paths.  He  assumed  two  simple  flow  geometries;  flow 
through  a  series  of  parallel  circular  ducts  and  flov;  through  a 
series  of  parallel  flat  ducts.  Biot  derived  equations  describing 
the  pore  fluid  flow  in  these  two  simplified  tore-  geometries. 


Flow  conditions  in  the  flat  duct  are  pictured 

schematically  in  Figure  2.1.  The  duct  is  assumed  to  have 

infinite  width  and  a  height  2a.  Flow  is  considered  over  a  unit 

width  of  duct  having  a  cross  sectional  area  of  2a.  Biot  (1956) 

expresses  the  ratio  of  viscous  shear  force  2t,  (where  t  is  Lhe 

viscous  shear  stress)  to  the  average  fluid  velocity  relative  to 

the  duct  walls,  W  .,  as 

ai 


(2-11) 
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The  average  relative  velocity  is  given  by 


W  .  =  U  .  -  u.  (2-12) 

ai  ai  1 

where  Ua^  is  the  average  absolute  fluid  velocity  in  the  duct  and 
u1  the  velocity  of  the  duct  walls,  y  is  the  dynamic  viscosity 
(dyr.e-sec/cm^ )  and  F(<)  is  a  complex  function  of  frequency  which 
describes  the  increase  in  viscous  shear  with  frequency.  Biot 
(1956)  gives  an  equation  for  F  of 


F ( z )  =  i  2  tanh(z) - 

1  -  -  tanh(2) 


(  2-13) 


where 


2=1 


1/2 


(2-14) 


and 


(2-15) 


Thus#  <  is  a  function  of  the  duct  dimension,  a,  the  dynamic 
viscosity,  y,  the  fluid  density,  p^,  and  the  excitation 
frequency  u .  The  friction  force  per  unit  volume  of  pore  fluid, 
from  liquation  2-8,  is  given  by 


D. 

l 


2t 

2a 


(  2-16) 


Substitution  of  Equation  2-11  into  2-in  gives 


(2-17) 


for  the  flat  duct. 

Flow  conditions  in  the  circular  duct  are  pictured 

schematically  in  figure  2.2.  In  this  case,  a  represents  the 

2 

radius  of  the  duct  having  a  cross  sectional  area  of  ^a  .  Biot 
(  1956)  gives  the  ratio  of  viscous  shear  force,  27rax,  to  the 
average  relative  fluid  velocity,  W  •,  as 


2iraT 

W  . 
ai 


8nyF  ) 


For  the  circular  duct, 


(2-18) 


F  (tc )  =  i  - -  (2-19) 

1  - 

where  k  is  given  in  Equation  2-15,  with,  a,  now  equal  to  the  radius 
of  the  circular  duct  and  the  function  T(ic)  given  by 


T(k) 


J0  UV2  K) 


( 2-20) 


Jc  represents  Bessel's  function  for  n  =  0  and  is  given  by 


Jo(x) 


6 

x 


2304 


(  2-21) 
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where 


x 


tc 


(  2-22) 


The  friction  force  per  unit  volume  of  pore  fluid  in  the  circular 
duct  D ^  from  Equation  2-8,  is  given  by 


D. 


2-rraT 


ira 


(2-23) 


Substitution  of  Equation  2-18  into  2-23  gives 


(2-24) 


for  the  circular  duct. 


Biot's  Approximations  for  Viscous  Shear 

Biot  (1962b)  derived  an  approximation  for  the  viscous 
friction  correction  factor,  F(k),  for  the  case  of  a  flat  duct, 
over  a  limited  range  of  frequencies.  An  alternative  derivation  to 
that  used  by  Biot  is  presented  h^re.  For  the  flat  duct,  the 
viscous  correction  factor  is  giv  .  by  Equation  2-13.  If  tanh(z) 
is  expressed  as 

13  2  5 

tanh(z)  =z-fZ+j^z~  ....  (2-25) 
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then 


z  tanh(z)  *=  7}  -  ^  z4  +  z3  -  . .  . .  (2-26) 

and 

1  -  tanh(z)  =  ^  z^  -  33"  2 ^  +  .  (2-27) 

Substitution  of  Equations  2-14,  2-26  and  2-27  into  Equation  2-13 
gives  an  approximation  for  the  flat  duct  of 


( 2-28) 


when  higher  order  terms  above  are  dropped.  This  agrees  with 
Biot's  approximation.  Biot  indicates  that  the  above  approximation 
is  accurate  to  within  5%  for  values  of  frequency  of 


0)  < 


(  2-29) 


and  represents  a  reasonable  approximation  of  F(k)  up  to  a 
frequency  of  about 


<x  < 
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a  P. 


(  2-30) 


2 

For  a  pore  size  of  a  :  .  005  cm  and  water  at  15°C  (p/p^  =  .013  cm  /sec) 
Biot  indicates  that  Equation  2-27  is  accurate  up  to  approximately 
300  Hz.  Mengi  and  McNiven  (1977)  determined  that  Biot's 
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approximation  given  in  Equation  2-27  is  accurate  even  for  the 
abrupt  stress  changes  at  a  wavefront. 

A  similar  approach  can  be  used  to  obtain  an  approximation 
of  the  viscous  correction  factor  for  the  circular  duct  given  by 
Equation  2-19.  Evaluating  T(k)  according  to  Equations  2-20  and 
2-21  gives 


TIO  =  j  i«  +  K*3  '  5f  i*5  +  ••••  l2-31> 

Substitution  of  Equation  2-31  into  2-19  gives  the  approximation 
for  F(k)  in  the  case  of  the  circular  duct  of 


F(k)  rx  1  +  (ix^)  +  .... 


(2-32) 


when  higher  order  terms  above  k  are  dropped. 


Biot  (1956)  indicates  that  there  is  an  approximate 
equivalence  between  the  circular  duct  and  the  flat  duct,  when  the 
height  of  the  flat  duct  equals  3/4  of  the  diameter  of  the 
circular  duct.  Applying  this  equivalence  to  Equations  2-29  and 
2-30  indicates  that  Equation  2-32  for  the  circular  duct  is 
reasonably  accurate  up  to  frequencies  of  nearly  twice  that  of  a 
flat  duct  if  the  respective  diameters  and  heights  are  equal. 

The  expressions  for  the  viscous  friction  correction 
factors  from  Equations  2-28  and  2-32  for  the  flat  and  circular 
ducts  respectively  can  be  substituted  into  Equations  2-17  and 
2-24  to  give  the  friction  force  per  unit  volume  of  pore  fluid  in 
each  type  of  duct.  For  the  flat  duct  this  yields 

D  =  w.  +  xk2  W.  (2-33) 

1  a  1  5a  1 
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and  for  the  circular  duct 


D. 

i 


(2-34) 


As  discussed  previously,  the  actual  relative  velocity  between  the 
skeleton  and  the  pore  fluid,  W^,  is  given  by 


W.  =  -  w. 

i  n  l 


( 2-35) 


Substitution  of  Equations  2-35  and  2-15  into  Equation  2-33  gives 
the  friction  force  per  unit  volume  of  pore  fluid  for  the  flat  duct 
as 


Di  "i  +  5n  iwWi  (2-36) 

na 

for  an  input  excitation  of  frequency  u> ,  the  excitation  velocity 
can  be  expressed  as 


w 


i 


( 2-37) 


where  w^  is  the  amplitude  of  the  excitation  velocity. 
Differentiating  Equation  2-37  with  respect  to  time  gives  the 
apparent  relative  acceleration,  w  as 

"  .  —  i  out  .  • 

wi  =  la)Wi  e  =  (2-38) 

Thus,  Equation  2-36  for  the  flat  duct  can  be  rewritten  as 
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E  2-39 ) 


For  the  circular  duct  substitution  of  Equations  2-35  and 
2-15  into  Equation  2-34  gives 


»i  *  ^ 


na 


pf  .  • 

+  —  lww . 

3n  i 


(  2-40 


Finally,  substitution  of  Equation  2-38  into  2-40  gives 


D  5=  _Ji_  w  + 

1  2  1 

pf  .. 

3n  Wi 

na 

Even  though  Equations  2-39  and  2-41  were  derived  assuming 
a  single  harmonic  forcing  function,  it  can  be  demonstrated  that 
they  should  apply  for  any  arbitrary  excitation.  Since  any 
arbitrary  excitation  function  can  be  decomposed  into  an  infinite 
series  of  harmonic  functions,  each  of  these  excitations  can  be 
treated  by  the  above  equations  and  the  total  friction  force 
obtained  by  superimposing  the  contribution  from  each  harmonic 
excitation  of  the  decomposed  total  input. 

Application  of  Darcy's  Law 

Biot's  equations,  based  on  thermodynamic  theory,  can  be 
related  to  Darcy's  flow  Law.  Darcy's  Law  assumes  that  the  flow 
velocity  is  proportional  to  the  hydraulic  gradient.  For  an 
arbitrary  gradient,  the  equivalent  permeability  coefficients  can 
be  derived  analytically  for  both  the  flat  and  circular  ducts. 
These  permeability  coefficients  can  then  be  compared  to  Biot's 
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flow  equations.  The  advantage  of  incorporating  Darcy's  Law  into 
the  Biot  equations  is  that  it  permits  use  of  parameters  which  are 
readily  obtained  from  laboratory  tests. 

The  flow  through  a  flat  duct  of  length  1  under  the 
pressure  gradient  p^  -  P2  is  pictured  schematically  in  Figure 
2.3.  The  heignt  of  the  duct  is  2a  and  a  unit  width  of  flow  is 
considered.  At  an  arbitrary  distance,  r,  from  the  center  of  the 
duct,  force  equilibrium  conditions  can  be  established  as 

Yfh1(2r)  -  Y f h 2  { 2 r )  +  y  §7  *21  =0  (2_42) 


where  y ^  is  the  unit  weight  of  fluid,  h^  and  h2  are  the  pressure 
heads  at  either  end  of  the  duct,  u  is  the  dynamic  viscosity  and  v 
is  the  fluid  velocity  at  distance  r.  The  first  two  terms  of 
Equations  2-42  represent  the  pressure  forces  P-^  and  P2  acting  on 
the  area  2r  at  the  left  and  right  ends  of  the  duct  respect ivev 
The  last  term  is  the  viscous  friction  force  acting  along  the  tw 
surfaces  at  a  distance  r  from  the  centerline  of  the  duct. 
Rearranging  Equation  2-42  as 

Yf  (h.  -  h,)  =  -y  P-  £  (2-43) 

'  f  1  2  dr  r 

and  substitution  of  the  hydraulic  gradient  i  given  by 


i  = 


(2-44) 


into  Equation  2-43  gives 
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dv  = 


P 


rdr 


(2-45) 


integrating  both  sides  of  Equation  2-45  from  o  to  r  gives  the 
fluid  velocity  at  the  distance  r  as 


v  = 


2 

r 


+  C 


(2-46) 


At  the  duct  wall  the  fluid  velocity  must  drop  to  zero  and 


C  = 


(2-47) 


Combining  Equations  2-46  and  2-47  yields 


v 


(2-48) 


The  average  velocity  in  the  duct,  V  #  is  defined  as 


a 


o 


where  Q  is  the  total  flow  through  the  duct  of  area  A  (equal  to 
2a).  Substitution  of  Equation  2-48  into  the  integral  of  2-49  and 
manipulation  gives 


v 


i 


(2-50) 


This  is  the  actual  average  fluid  velocity  relative  to  the  duct. 
The  duct  represents  only  the  pore  volume  of  the  soil.  As 
explained  in  conjunction  with  Equation  2-5,  the  actual  discharge 
velocity  w^ ,  or  apparent  relative  fluid  velocity,  is  given  by 


w . 

i 


nv  = 
a 


na  y j 
3y 


(2-51) 


Darcy's  Law  states  that 


w.  =  ki 

i 


(2-52) 


where  k  _ _  the  coefficient  of  permeability  (cm/sec).  From 
Equations  2-51  and  2-52  the  coefficient  of  permeability  for  the 
flat  duct  can  be  expressed  as 


(2-53) 


A  similar  procedure  can  be  followed  for  obtaining  the 
average  apparent  flow  velocity  through  a  circular  duct.  This 
solution,  known  as  the  Hagen-Poiseu 1 le  flow  law  is  given  by 


2 

na  y 


w .  = 
i 


8y 


1  i 


(2-54) 


where  a  is  the  radius  of  the  circular  duct  (e.g.  Bowles,  1979). 
From  Equations  2-52  and  2-53  the  coefficient  of  permeability  fo 
the  circular  duct  is  given  as 


(2-55) 


Biot's  equations  for  viscous  friction  (2-39  and  2-41)  can  now  be 
expressed  in  terms  of  the  coefficient  of  permeability  given  in 
Equations  2-53  and  2-55.  For  the  flat  duct, 


Tf  •  pf  •• 

D .  =  -s—  w  .  +  c —  w . 
l  k  i  5n  l 


(2-56) 


and  for  the  circular  duct 


Yf  .  Pf  .. 

Di  “  TT  wi  +  3n  wi 


(2-57) 


Note  that  in  both  of  the  above  equations  when  the  input 
excitation  frequency  drops  to  zero  the  last  term  drops  out  and 
Biot's  equations  agree  exactly  with  Darcy’s  Law. 


Explicit  Form  of  Biot's  Equations 

The  explicit  friction  force  terms  from  Equations  2-56  and 
2-57  for  the  flat  and  circular  ducts  can  be  substituted  into 
Equation  2-10  for  the  pore  pressure  gradient  giving 


f  i 

—  (1  +  5>  w.  +  pfu. 


( 2-58 ) 


for  the  flat  duct,  and 


( 2-59) 

for  the  circular  duct. 

The  inclusion  of  the  explicit  approximations  for  viscous 
friction  in  Equations  2-58  and  2-59  is  in  the  form  of  a  viscous 
shear  given  by  Darcy's  Law  in  the  last  term  plus  an  added 
increment  of  inertial  resistance  associated  with  the  apparent 
acceleration  of  the  pore  fluid  relative  to  the  skeleton.  For  the 
flat  duct  this  increment  is  1/5  and  for  the  circular  duct  it  is 
1/3.  The  increment  depends  only  on  the  pore  shape  and  not  the 
size. 

As  discussed  previously  in  the  subsection  on  Biot's 
equations,  it  was  explained  that  the  assumed  flow  paths  in  the 
form  of  flat  and  circular  ducts  were  a  poor  approximation  of  the 
actual  paths  in  saturated  soil.  In  order  to  study  the  influence 
of  more  restrictive  flow  paths,  Equations  2-58  and  2-59  can  be 
generalized  in  the  following  form 

(2-60) 

where  r  represents  the  additional  inertia  or  mass  increment 
factor  which  is  a  function  of  the  flow  path  geometry. 

TP DAP  FORMULATIONS 

The  results  of  the  previous  analytic  development  are 
incorporated  into  the  two  phase  finite  element  code  TPDAP.  Based 
on  the  governing  differential  equation  for  the  bulk  mixture 
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(Equation  2-7)  and  for  the  pore  fluid  (Equation  2-60),  global 
equilibrium  equations  for  two  phase  media  are  established.  These 
global  equilibrium  equations  are  presented  in  discretized  form  in 
both  the  space  and  time  domain.  TPDAP  uses  either  Newmark's  6 
method  or  Wilson's  8  method  for  time  integration. 


Notation 

Positive  signs  are  used  for  elongation  and  tension.  A 
comma  denotes  differentiation  with  respect  to  the  subsequent 
indices  and  a  superposed  dot  denotes  time  rate. 


total  stress 


effective  stress 


pore  fluid  pressure 

solid  skeleton  displacement 

apparent  fluid  displacement  relative  to  solid 
skeleton 

solid  skeleton  strain 

solid  skeleton  volumetric  strain 
fluid  phase  volumetric  strain 

nodal  solid  skeleton  displacement  vector  at  the 
element  degrees  of  freedom 


nodal  apparent  relative  fluid  displacement 
vector  at  the  element  degrees  of  freedom 
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{T} 
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TT 


[<3] 

[Dep] 
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m 


6,  Y 


e 


r 


s  applied  boundary  traction  itotal  stress 
:  specified  boundary  pore  fluid  pressure 
:  inverse  of  permeability  matrix 
:  elasto-plastic  stress-strain  natrix 
:  porosity 

:  bulk  modulus  of  pore  fluid 
:  bulk  modulus  of  solid  grain 
:  bulk  modulus  of  soil/water  mixture 
s  bulk  mass  density  of  mixture 
s  fluid  mass  density 
j  parameters  of  Newmark's  6  method 
:  parameter  in  Wilson's  6  method 

:  Kronecker's  delta 

f.  •  =  0  if  i  f  j 
13 

5 .  .  =  1  if  i  -  j 
ID 

:  mass  increment  factor 
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Field  Equal.  ions  i_n  Two-Phase  Dynami  cs 

Effective  stress  strain  relation*. 

tiy  virtue  of  Terzaghi's  effective  stress  law,  the  total 
stress  is  expressed  r.s  the  sum  of  the  effective  stress  and  the 
pore  water  pressure 


o 


ij 


=  r 


.  .  +  6  .  .  u 
ID  ID 


(2-61) 


Strain-displacement  relation: 

Under  the  assumption  of  small  displacement  theory, 


E 


+ 


(2-62) 


Continuity  equation: 

Based  on  the  conservation  of  mass,  the  coupled  continuity 
equation  of  flow,  as  derived  by  Kim  (1982),  is  given  by 


dit 


n)  dc  +  n  dt. 


K 

m 


(2-63) 


where 


K 


m 


K  K 


w 


K 

g 


kT 

w 


(  2-64  ) 
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Solid  skeleton  constitutive  law: 


neP 


i  jkLdek)t 


(  2-65) 


Equation  of  motion  for  the  bulk  mixture: 

The  differential  equation  of  motion  governing  the  bulk 
mixture  was  derived  earlier  in  this  section  and  is  expressed  as 
follows : 


o..  .  =  pu.  +  prw. 

l] , D  H  1  l 


(2-7) 


Equation  of  motion  for  the  pore  fluid: 


The  generalized  form  of  Biot’s  equation  was  given  by: 


(1  +  r)  w.  +  P  u.  + 
n  ill 


f  • 

— T—  W  . 

k  l 


2-60) 


Equation  2-60  was  derived  under  the  assumption  of  isotropic 
permeability.  In  an  anisotropical  ly  permeable  porous  medium,  the 
coefficient  of  isotropic  permeability,  k,  can  be  replaced  by  the 
symmetric  anisotropic  permeability  matrix,  [k],  given  by 
Scheidegger  (1957)  as 


XX 

xy 

xz 

k 

k 

xy 

yy 

yz 

k 

k 

xz 

yz 

zz 

2  4 


[k]  = 


(2-66) 


and  the  generalized  torra  of  Clot's  equation  in  anisotropieally 
permeable  medium  can  be  expressed  as 


n  '3- 


+  r>  "i  +  °£ui  +  wi 


( 2— 67  > 


where 


[q]  = 


(2-68) 


Discret  ized  Global  Equilibrium  Equation 

Two  global  equilibrium  equations  are  derived,  first  in 
terms  of  field  variables  and  then  discretized  by  the  nodal  variables. 


The  first  equation  equates  the  total  internal  stresses 
plus  the  inertia  forces  to  the  applied  boundary  tractions. 
Letting  the  solid  skeleton  movement  be  the  virtual  displacement, 
6u,  the  following  global  equilibrium  equation  for  the  bulk 
mixture  is  established, 


L 


■  Ij 


(6e)  (o}dv  =  /{6u}1{T}ds 


•  L 


{Su}1  p{(i}dv 


-L 


(6u)  1  pf  {tt}  dv 


(2-69! 
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where  6e  is  the  virtual  strain  corresponding  to  the  virtual 
displacement  6u. 


The  second  equation  equates  the  applied  pore  pressure  on  the 
boundary  to  the  internal  pore  pressure  plus  the  flow  resistance 
force  plus  the  inertia  force  on  the  pore  fluid.  Taking  the 
apparent  relative  fluid  movement  as  the  virtual  displacement,  <5w, 
the  internal  virtual  work  done  by  the  pore  pressure  should  be 
equal  to  the  external  virtual  work.  That  is, 


/  (6W.  .)Tirdv  =  /  {6w}T  nds  -  /  {6w}T[q]{w) 

1 ' 1  Jv  -'v 


dv 


-  f  (5w}T  pf{ii}dv  -  f  {( Sw}T  ^  pf(w)dv  (2-70) 

Jsj  J v 


Field  variables  can  be  discretized  into  nodal  values  withir.  the 
finite  element. 


(u) 

U) 

(w) 


w . 

1,1 


=  [N]  { ule 
=  [B]  { u)e 
=  [N]  {w)e 
=  {1}T[b]  {w)e 


(2-71  ) 


where  {1}=  <111000> 


Replacing  the  field  variables  in  Equations  2-69  and  2-70 
by  the  discretized  nodal  variables  using  Equation  2-71  gives  the 
following  global  equilibrium  equations  at  time  step  n; 
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or 


MU  +  D  U  +  K  AU  =  P  -  R  , 
n  n  n  n-J. 


( 2-72) 


where 
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m 
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K, 


n-1 


n-1 


=  I  J [N]T  p[N]dv 

-  £  J  [N]T  pf  [N]dv 

=  I  J [B]T  [DeP]  [B]dv 

=  l  f [B]T  {1)K  {l}T[B]dv 
Jw 

=  I  j [B]T{o'n_1}dv 
=  Z  f  [B]T  {  1) TT  , dv 
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Now  introducing  Newmark's  B  method  and  Wilson's  9 
method  for  time  integration,  fc’quation  2-72  can  be  expressed  in 
the  following  form; 


K 


t 


n 


(2-73) 


where 


K  + 


M  + 


6t‘ 


(2-74) 
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SECTION  3 


TWO  PHASE  DYNAMIC  ANALYSIS  - 


SATURATED  POROUS  ELASTIC  SKELETON 


INTRODUCTION 

The  next  two  sections  present  a  series  of  parameter 
studies  in  which  saturated  two  phase  porous  materials  are 
subjected  to  dynamic  uniaxial  loadings.  These  parametric 
calculations,  using  the  TPDAP  code,  were  arranged  so  that  the 
influence  of  the  various  material  and  loading  parameters  on  the 
overall  response  could  be  isolated  and  analyzed  individually. 
Three  different  soil  skeletons  were  modeled.  In  this  section  a 
linear  elastic  skeleton  was  studied.  In  the  following  section  a 
bilinear  hysteretic  skeleton  was  modeled,  followed  by 
calculations  on  an  actual  sand  from  Enewetak  Atoll. 

In  this  section  three  different  loadings  are  applied;  a 
step  (or  Heaviside)  loading  pulse,  a  triangular  loading,  and  a 
fifth  order  simple  polynomial  loading.  All  loadings  applied  a 
peak  stress  of  5  ksi.  other  parameters  investigated  in  the 
elastic  calculations  include  a  comparison  between  two  phase  and 
total  stress  (one  phase)  analysis,  the  influence  of  mass 
discretization  techniques,  the  influence  of  Biot's  additional 
inertia  term,  and  the  influence  of  permeability,  numerical 
damping,  and  surface  drainage  conditions. 
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CALCULATIONAL  AND  MATERIAL  PARAMETERS 


Loading  Functions 

The  three  5  ksi  loading  time  histories  used  in  this 
section  are  shown  in  Figure  3.1.  The  rise  time  for  all  loadings 
was  0.1  msec.  The  triangular  and  polynomial  loadings  had  a 
positive  phase  duration  of  20.1  msec.  The  pressure  decay  of 
polynomial  loading  is  given  by 


P(t)  =  p. 


(t  -  tr) 


(3-1) 


where  pQ  is  the  peak  pressure,  t  the  time,  tr  the  rise  time,  and 
tQ  the  decay  time.  The  total  impulse  under  the  polynomial 
loading  is  one  third  of  the  impulse  under  the  triangular  loading. 


’■'in ite  Element  Parameters 

All  the  calculations  throughout  this  study  used  the 
finite  element  mesh  shown  in  Figure  3.2.  The  mesh  consists  of 
200  elements  with  a  total  depth  of  100  ft.  The  elements  increase 
in  thickness  with  increasing  depth,  starting  with  a  thickness  of 
0.2  ft  at  the  loading  surface  and  expanding  to  0.8  ft  at  the  100 
ft  depth.  A  constant  time  step  of  0.05  msec  was  used  in  all 
ca leu  La  t ions . 

A  number  of  parameters  were  standardized  for  the 
calculations  in  this  section.  These  were  held  constant  in  most 
of  the  calculations,  except  when  the  influence  of  varying  any  one 
of  these  parameters  was  examined.  The  most  suitable  value  of 
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each  parameter  was  not  known  beforehand;  thus,  the  sensitivity  to 
each  was  examined  as  part  of  this  parameter  study. 

In  most  calculations,  the  loading  pressure  function  was 
partitioned  on  the  ground  surface  in  proportion  to  the  areal 
ratio  of  the  solid  and  fluid  constituents.  In  other  words,  n%  of 
the  total  force  on  the  loading  boundary  was  applied  to  the  pore 
water  and  the  remaining  (1  -  n)%  was  applied  to  the  soil 
skeleton.  For  the  porosity  of  35%  and  peak  stress  of 

5  ksi,  the  peak  stress  applied  to  the  pore  water  at  the  loading 
boundary  was  1750  psi  while  the  peak  effective  stress  applied  to 
the  soil  skeleton  was  3250  psi. 

In  the  discussion  of  Biot’s  equations  of  Section  2,  an 
additional  inertia  term  (Equation  2-60)  was  added  to  TP DAP  to 
account  for  the  increased  friction  forces  which  develop  at  higher 
frequencies.  In  most  of  the  calculations  in  this  section,  r  is 
given  a  value  of  zero;  i.e.  there  is  no  added  friction  force  due 
to  relative  acceleration  between  the  skeleton  and  pore  fluid. 

The  only  friction  force  between  the  two  is  proportional  to  the 
relative  velocity  and  is  given  by  Darcy's  law. 

For  high  frequency  loadings,  such  as  used  in  this  study, 
oftentimes  oscillations  develop  at  the  wavefront.  To  minimize 
the  potential  for  oscillations,  numerical  damping  was  introduced 
into  most  calculations  in  this  section.  In  this  study  Newmark's 
method  for  time  integration  was  selected,  with  y  =  1.2,  which 
results  in  relatively  heavy  damping.  The  corresponding  value  of 

6  -  0.7225  calculated  from 


6 


(Y  +  0.5)2 

4 


(  3-2) 


was  used. 


Archer  (1963)  demonstrated  that  the  conventional  lumped 
treatment  of  mass  in  finite  element  calculations  can  result  in 
large  errors  when  significant  inertial  forces  are  present.  More 
accurate  solutions  can  be  obtained  by  use  of  the  consistent  mass 
matrix  technique  described  by  Archer,  in  which  the  matrix  is 
consistent  with  the  actual  distribution  of  mass  in  the  medium. 

The  use  of  the  consistent  mass  technique  results  in  a  symmetric 
mass  matrix  rather  than  the  simple  diagonal  mass  matrix  of  the 
lumped  mass  approach.  In  most  of  the  calculations  in  this  study 
the  consistent  mass  technique  is  used. 

Materia  1  Parameters 

The  material  properties  used  in  the  linear  elastic 
analysis  are  summarized  in  Table  3.1.  The  skeleton  properties  are 
representative  of  a  uniform  medium  dense  uncemented  quartz  sand. 
The  elastic  moduli,  From  Blouin  and  Kim  (February,  1984)  are  the 
secant  moduli  at  1.5%  strain,  obtained  from  laboratory  data  on  5 
typical  sands.  The  undrained  bulk  modulus  of  the  soil  water 
mixture,  K^,  is  obtained  from  the  Wood  equation  and  represents 
the  modulus  of  the  soil  water  mixture  with  no  stiffness 
contribution  from  the  soil  skeleton  (for  further  discussion  see 
Blouin  and  Kim,  February  1984  and  Richart  et  al.,  1970).  Based 
on  the  Wood  equation,  a  compre ss iona 1  wavespeed  of  5161  ft/sec  is 
computed . 

The  simplest  model  for  the  composite  undrained 
compressibility  of  the  saturated  soil  is  the  decoupled  model 
(Blouin  and  Kim,  February,  1984).  The  decoupled  model  assumes 
that  the  stress  is  resisted  by  the  stiffness  of  the  soil  skeleton 
acting  in  parallel  with  the  stiffness  of  the  soil-water  mixture 
from  the  Wood  equation.  The  resultant  bulk  and  constrained 
moduli,  K^  and  are  simply  the  sum  of  the  mixture  modulus  and 

the  bulk  and  constrained  moduli  of  the  soil  skeleton 
respectively.  In  the  undrained  case,  the  TP DAP  continuity 


equation  combined  with  Tertaghi's  effective  stress  law  results  in 
the  decoupled  modulus  equations. 


Transmitting  Boundary 

In  order  to  prolong  the  calculation  time  available  with 
the  mesh  of  Figure  3.2,  a  transmitting  boundary  was  employed  on 
the  base  of  the  botcom  element.  In  establishing  this  boundary, 
it  was  assumed  that  there  was  no  relative  motion  between  the 
pore  fluid  and  solid  skeleton  at  the  boundary.  This  condition  is 
expressed  as 


w . 
i 


=  0 


(3-3) 


where  w^  is  the  apparent  fluid  displacement  relative  to  the  soil 
skeleton  at  the  boundary.  Thus,  the  saturated  soil  is  assumed  to 
have  no  drainage  at  the  boundary  and  the  equations  for  the 
undrained  decoupled  modulus  of  Table  3.1  apply.  This  implies 
that  beneath  the  boundary  one  phase  total  stress  analysis  is 
valid,  where  the  compress iona 1  wavespeed 


(3-4) 


from  Table  3.1  is  used. 

The  above  assumptions  permit  use  of  the  transmitting 
boundary  developed  by  Lysmer  and  Kuhlemeyer  (1969)  where  the 
stress  wave  energy  on  the  boundary  is  absorbed  by  dashpots  at  the 
element  nodes.  The  total  normal  stress  on  the  boundary,  o  ,  is 
given  by 
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where  un  is  the  absolute  velocity  of  the  skeleton  normal  to  the 
boundary.  In  order  for  the  energy  dissipated  in  the  dashpots  to 
just  equal  the  energy  in  the  compress ional  wave,  the  dashpot 
constants  are  proportional  to  the  undrained  impedence  pC^. 

Because  there  is  relative  motion  between  the  skeleton  and 
the  pore  fluid,  the  above  assumptions  do  not  entirely  eliminate 
reflections  frcm  the  boundary.  Development  of  truly 
non-reflecting  boundaries  for  two  phase  media  is  beyond  the  scope 
of  this  work.  However,  use  of  the  above  assumptions  in  the 
elastic  two  phase  calculations  greatly  reduced  reflections  from 
the  boundary. 


HEAVISIDE  LOADING 

The  Heaviside  loading  from  Figure  3.1  was  used  in  the 
initial  series  of  parametric  calculations  because  it  represented 
a  simple  loading  and  permitted  study  of  pore  pressure  dissipation 
toward  the  ground  surface. 


One  Phase  Calculations 


Two  conventional  one  phase  calculations  were  performed 
under  the  Heaviside  loading  with  mass  discretization  using  the 
consistent  and  lumped  techniques.  These  are  compared  in  Figure 
3.3  in  the  form  of  pore  pressure  profiles  at  times  of  2,  6,  and 
10  msec.  The  material  properties  from  Table  3.1  were  used,  with 
the  moduli  given  oy  the  undrained  decoupled  moduli,  K ^  and 
The  results  of  the  one  phase  total  stress  calculations  were 
converted  to  pore  pressure  for  comparison  to  the  pore  pressure 
profiles  obtained  in  the  later  two  phase  calculations. 
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The  total  vertical  stress/  a  ,  in  uniaxial  loading  is  given 
by 


a 


V 


M,e 

d  v 


(3-6) 


where  ev  is  the  vertical  strain  and  is  the  decoupled 
constrained  modulus  from  Table  3.1.  As  described  by  Blouin  and 
Kim  (February,  1984),  the  effective  stress,  a  v' ,  in  uniaxial 
loading  is  given  by 
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(  3-7  ) 


where  Ms  is  the  constrained  modulus  of  the  skeleton.  According 
to  Terxaghi's  effective  stress  law,  the  total  stress  is  comprised 
of  the  effective  stress  plus  the  pore  water  pressure,  n, 
as 


a  =  n  +  o  ' 
v  v 


(3-8) 


The  ratio  of  the  pore  pressure  to  the  total  vertical  stress  in  an 
undrained  uniaxial  loading  is  obtained  by  solving  Equation  3-8 
for  it  and  substituting  the  values  of  total  and  effective  stress 
from  Equations  3-6  and  3-7  into  the  ratio  to  give 


7T 


(  3-9) 
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Using  the  values  of  and  Mg  from  Table  3.1  gives  the  ratio  as 


—  =  0.892  ( 3-10) 

°v 

Thus,  the  pore  pressure  profiles  of  Figure  3.3  were  obtained  by 
multiplying  the  total  vertical  stress  from  the  one  phase 
calculations  by  the  ratio  of  pore  pressure  to  total  stress  given 
in  Equation  3-10. 

As  demonstrated  in  Figure  3.3,  there  is  no  significant 
difference  between  the  one  phase  calculations  using  the  lumped 
mass  technique  and  those  using  the  consistent  mass  technique. 
Because  the  mesh  used  in  these  calculations  is  so  fine,  the 
errors  introduced  by  the  lump  mass  technique  into  the  one  phase 
calculation  are  negligible. 

TWO  PHASE  CALCULATIONS 

Lumped  Mass  vs .  Consistent  Mass  Calculations 

The  initial  set  of  two  phase  calculations  under  the 
Heaviside  loading  compared  the  computed  response  using  the  lumped 
mass  technique  with  that  of  the  consistent  mass  technique.  In 
contrast  to  the  one  phase  results,  there  was  a  significant 
difference  between  the  two  techniques  in  the  two  phase 
calculations.  This  difference  is  illustrated  in  the  stress 
profile  comparisons  at  2,  5,  10  and  14  msec  shown  in  Figure  3.4. 

As  discussed  previously,  the  5  ksi  Heaviside  loading  is 
applied  at  the  surface  and  is  partitioned  between  the  soil 
skeleton  and  the  pore  water  at  the  loading  boundary.  For  the 
assumed  porosity,  35%  or  1750  ps i  is  applied  to  the  pore  water 
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and  the  remaining  65%  or  3250  psi  is  applied  to  the  solid 
skeleton.  Despite  this  partitioning  at  the  surface,  Figure  3.4 
shows  that  within  several  feet  from  the  surface  the  load  is 
redistributed  between  the  two  phases  so  that  the  pore  water 
carries  about  90%  of  the  total  load,  as  specified  by  Equation 
3-10  for  the  undrained  condition.  Because  of  the  fluid  friction 
between  the  pore  water  and  the  skeleton,  the  stress  distribution 
applied  at  the  surface  is  very  rapidly  transformed  to  that  which 
would  result  from  an  undrained  loading. 

The  permeability  of  the  soil  skeleton  permits  flow  of 
pore  water  toward  the  surface.  This  flow  results  from  the  pore 
pressure  gradient  between  the  high  pressure  in  the  undrained 
region  behind  the  wavefront  and  the  1750  psi  pore  pressure 
applied  at  the  ground  surface.  The  pore  water  migration  from  the 
high  pressure  region  toward  the  low  pressure  boundary  region 
causes  dissipation  of  the  pore  pressure.  As  the  pore  pressure 
decreases,  the  intergranular  or  effective  stress  increases  as  the 
skeleton  compresses  and  assumes  an  increasing  portion  of  the 
total  stress.  At  all  depths,  the  sum  of  the  pore  pressure  and 
the  effective  stress  equals  the  total  stress.  The  pore  pressure 
dissipation  advances  downward  from  the  surface  with  increasing 
time.  The  rate  of  advancement  is  a  function  of  the  permeability 
and  the  pore  pressure  gradient  as  in  the  usual  consolidation 
process.  At  14  msec,  the  pore  pressure  dissipation  front  has 
reached  a  depth  of  about  15  ft. 

At  the  wavefront  in  Fiure  3.4,  the  profile  computed  using 
the  lumped  mass  matrix  is  considerably  more  smeared  than  that 
computed  from  the  consistent  mass  matrix.  In  Equation  2-72  the 
first  matrix  multiplication  term 
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expresses  the  inertial  resistance  of  the  bulk  mixture  and  the 
pore  water.  The  inertial  force  vector  for  the  soil  water  mixture 
(Im}  is  given  by 


' V  -  fU'V  *  Cmc]  (3-ni 

and  the  inertial  force  vector  for  the  pore  fluid  {1^}  is  given  by 

(If)  =  [Mc]T  {un>  +  [Mf]{wnl  (3-12) 

The  consistent  mass  matrix  includes  all  the  mass  (M)  terms  in  the 
above  matrix.  However,  the  lumped  mass  matrix  includes  only  the 
Mm  and  terms  in  the  above  matrix.  In  addition,  the  lumped 
mass  matrices  and  are  diagonalized,  which  introduces  an 
additional  error.  The  effect  of  neglecting  the  coupling  mass 
matrix  term,  Mc,  drops  the  second  inertial  force  term  out  of 
Equation  3-11  and  the  first  inertial  force  term  out  of  Equation 
3-12.  In  most  problems  the  acceleration  of  the  skeleton,  un,  is 
much  greater  than  the  relative  acceleration  between  the  pore 
fluid  and  the  skeleton,  wn.  In  Equat.:on  3-11,  the  coupling  mass 
matrix,  Mc,  is  also  smaller  than  the  lass  matrix  for  the  mixture, 
Thus,  the  second  inertial  force  term  is  of  relatively  minor 
importance.  However,  in  Equation  3-12  the  coupling  mass  matrix 
is  associated  with  the  relatively  large  skeleton  acceleration;  and 
the  contribution  of  this  term  is  substantial. 

The  first  inertial  force  term  in  Equation  3-12 
corresponds  to  the  second  inertial  force  term,  PfU^,  in  the 
governing  differential  equation  of  motion  of  the  pore  fluid. 
Equation  2-60.  Inspection  of  Equation  2-60  shows  that  neglecting 
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this  term  will  significantly  increase  the  apparent  relative 
velocity  and  acceleration/  and  w^,  for  a  given  pore  pressure 
gradient.  This  increased  relative  motion  will  result  in 
unrealistically  high  energy  dissipation  and  smearing  of  the 
wavefront  as  shown  in  Figure  3.4.  In  general,  the  lumped  mass 
technique  is  inappropriate  for  use  in  two  phase  dynamic  analysis. 

Comparison  of  One  and  Two  Phase  Ca leu lat ions 

The  pore  pressure  profiles  from  the  one  phase  calculation 
using  the  consistent  mass  matrix  are  compared  to  those  from  the 
two  phase  calculation  in  Figure  3.5.  The  pore  pressure 
dissipation  near  the  surface  in  the  two  phase  calculation  is 
clearly  evident.  At  the  wavefront  the  two  phase  calculation  is 
smeared  considerably;  however,  this  is  not  a  result  of  numerical 
smearing  but  is  due  to  the  dissipation  of  energy  near  the  front 
resulting  from  the  relative  motions  induced  between  the  pore 
water  and  soil  skeleton.  The  effect  is  similar  to  that  caused  by 
numerical  damping,  but  in  this  case  is  a  real  physical  effect. 

Influence  of  Mass  Increment  Factor 

The  final  parameter  examined  in  the  Heaviside  loading 
calculations  was  the  influence  of  the  mass  increment  factor,  r, 
included  in  Equation  2-60.  As  explained  in  Section  2,  r  is  a 
factor  which  attempts  to  account  for  frictional  resistance 
between  the  pore  fluid  and  skeleton  in  excess  of  the  resistance 
accounted  for  by  Darcy's  law.  For  Biot's  evaluation  of  the 
circular  and  flat  ducts,  the  factor  r  depends  on  the  shape  of 
the  duct  but  is  independent  of  the  size.  The  value  of  r  for  the 
circular  duct  is  1/3  and  the  value  for  the  flat  duct  is  1/5. 

Figure  3.6  compares  pore  pressure  and  effective  stress 
profiles  for  two  phase  calculations  with  r  equal  to  zero  and 
r  equal  to  1/3.  For  an  r  of  1/3  there  is  a  slight  reduction  in 
the  smearing  at  the  wavefront.  The  additional  frictional 


43 


resistance  due  to  the  r  factor  results  in  somewhat  reduced 
relative  motion  and  energy  dissipation;  a  trend  toward  less 
smearing  of  the  wavefront  as  seen  in  the  undrained  one  phase 
calculation.  As  noted  in  the  discussion  of  Section  2,  the  actual 
pore  shapes  are  far  more  complicated  than  the  simple  models  used 
by  Biot,  so  realistic  values  of  r  raa'f  he  much  greater  than  those 
he  derived.  This  is  an  area  requiring  further  research. 


TRIANGULAR  LOADING 

One  Phase  vs.  Two  Phase  Calculations 


The  triangular  loading  function  of  Figure  3.1  was  used  in 
a  second  series  of  parametric  calculations.  With  a  rise  tine  of 
G. 1  msec  and  a  positive  phase  duration  of  20.1  msec,  this 
pressure  function  more  nearly  approximates  certain  airblast 
loadings.  The  first  series  of  calculations  using  the  triangular 
loading  pulse  compared  the  pore  pressure  response  from  one  and 
two  phase  calculations  as  shown  in  Figure  3.7.  In  the  one  phase 
calculation,  plotted  as  a  dashed  line,  the  pore  pressure  rises 
rapidly  to  a  peak  which  is  followed  by  a  linear  decay  which 
mirrors  the  linear  decay  of  the  loading  waveform.  As  explained 
in  conjunction  with  the  one  phase  calculation  for  the  Heaviside 
loading,  the  pore  pressure  is  obtained  by  multiplying  the  total 
stress  by  0.892  yiven  by  Equation  3-9.  The  solid  line  shows  the 
profiles  from  the  two  phase  calculation.  At  the  wavefront,  there 
is  significant  smearing  compared  to  the  one  phase  calculation 
and  some  attenuation  of  the  peak  stress.  Both  the  smearing  and 
additional  peak  stress  attenuation  are  primarily  due  to  the 
energy  dissipation  associated  with  the  relative  motion  between  the 
pore  fluid  and  the  skeleton.  Note  that  there  is  also  some 
attenuation  of  the  peak  stress  in  the  one  phase  profile  as  the 
wave  propagates  downward.  At  the  surface  the  applied  peak  stress 
is  4460  psi.  At  14  msec  and  a  depth  of  68  ft,  the  peak  stress 
has  attenuated  by  9%  to  about  4060  psi.  This  stress  attenuation 
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in  the  one  phase  calculation  is  due  to  the  numerical  damping.  The 
additional  peak  stress  attenuation  in  the  two  phase  calculation 
is  due  to  dissipation  by  fluid  friction.  At  14  msec  in  the  two 
phase  calculation  the  peak  stress  has  attenuated  by  an  additional 
9%  to  a  value  of  3660  psi  at  a  depth  of  64  ft. 

In  the  near  surface  region  the  two  phase  calculation 
shows  a  decrease  in  pore  pressure  resulting  from  partial  drainage 
toward  the  surface.  At  14  msecs  the  pore  pressure  dissipation 
front  has  reached  a  depth  of  13  ft,  slightly  less  than  that 
observed  for  the  Heaviside  loading.  This  slower  dissipation 
front  is  due  to  the  smaller  pore  pressure  gradients  under  the 
triangular  loading. 

Influence  of  Permeabi 1 i ty 

The  second  parameter  examined  in  the  series  of  triangular 
loading  calculations  was  the  influence  of  permeability  on  the 
material  response.  Figure  3.8  compares  the  pore  pressure  and 
effective  stress  profiles  at  various  times  between  a  soil  having 
the  standard  0.1  in/sec  permeability  and  a  soil  having  a 
permeability  of  0.001  in/sec.  The  higher  permeability  is  typical 
of  a  coarse  uniform  sand  of  high  permeability  while  the  lower 
value  is  typical  of  well  graded  sands  of  medium  permeability. 

There  are  significant  differences  between  the  two  calculations 
both  at  the  wavefront  and  in  the  near  surface  region.  As 
expected,  pore  pressure  dissipation  toward  the  surface  is 
severely  curtailed  in  the  lower  permeability  material.  At  14  msec 
the  pore  pressure  dissipation  front  has  reached  a  depth  of  less 
than  2  ft,  compared  to  about  13  ft  in  the  more  permeable  soil. 

The  wavefront  in  the  lower  permeability  soil  is 
significantly  less  smeared  than  that  in  the  higher  permeability 
soil.  The  lower  permeability  inhibits  relative  motion  between 
the  pore  fluid  and  soil  skeleton  at  the  wavefront,  thus, 
lessening  the  energy  dissipation  effects  due  to  fluid  friction. 
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Comparing  the  profiles  for  the  two  phase  calculation  in  the  low 
permeability  soil  to  the  one  phase  calculation  plotted  in  Figure 
3.7  shows  that  the  two  phase  calculation  is  nearly  identical  to 
the  one  phase  undrained  case.  Under  the  one  dimensional  linear 
elastic  assumptions  of  this  set  of  calculations,  the  one  phase 
calculation  adequately  models  the  response  of  this  medium 
permeability  soil  to  the  dynamic  loading. 


Inf luence  of  Mass  Increment  Factor 

The  third  parameter  investigated  in  this  series  of 
calculations  was  the  influence  of  the  mass  increment  factor  r. 
Under  the  Heaviside  loading,  use  of  an  r  value  of  1/3, 
corresponding  to  Biot's  circular  pores,  had  only  a  small 
influence  on  the  response  at  the  wavefront  (see  Figure  3.6). 
Figure  3.9  compares  the  effective  stress  and  pore  pressure 
profiles  for  calculations  using  r  values  of  zero  and  one.  An  r 
value  of  one  was  picked  to  accentuate  the  influence  of  r.  There 
is  less  smearing  of  the  wavefront  for  the  calculation  using  an  r 
of  one,  the  difference  being  considerably  greater  than  in  the 
previous  comparison  using  an  r  value  of  1/3.  As  was  mentioned 
previously,  increasing  the  r  value  results  ir.  effectively 
increasing  inertial  resistance  in  the  pore  fluid  while  reducing 
the  relative  fluid  motions.  However,  it  is  not  known  whether  an 
r  of  one  is  a  meaningful  value  for  this  parameter. 


Inf  luence  oj^  Damping 

The  next  parameter  studied  in  the  triangular  loading 
series  was  the  influence  of  the  Y  damping  on  the  response  of  the 
two  phase  material.  A  rather  heavy  Y  damping  of  1.2  was  used  as 
the  standard  throughout  most  of  this  study.  Figure  3.10  compares 
the  pore  pressure  and  effective  stress  profiles  calculated  using 
the  standard  damping  and  no  y  damping,  i.e.  y  -  0.5.  Behind  the 
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peak  stresses,  the  two  solutions  are  identical.  At  the 
wavefront,  the  undamped  profile  exhibits  a  sharp  jump  in  pote 
pressure,  as  compared  to  the  more  smeared  increase  in  pressure  ir. 
the  damped  case.  At  early  time,  the  undamped  stresses  exhibit 
oscillations  just  behind  the  wavefront.  Beyond  10  msec  (below  a 
depth  of  60  ft),  these  oscillations  no  longer  occur.  Finally, 
the  smearing  around  the  peak  stress  in  the  damped  calculation  is 
duplicated  in  the  undamped  calculation,  indicating  that  it  is 
indeed  caused  by  the  dissipation  resulting  from  relative  motion 
between  the  two  phases  near  the  wavefront  and  is  not  associated 
with  numerical  damping. 

Another  comparison  of  the  influence  of  damping  is 
illustrated  in  Figure  3.11  where  the  standard  damping  with  y  = 

1.2  is  compared  to  the  response  for  a  damping  of  0.85,  midway 
between  the  previous  comparisons.  There  is  little  difference 
between  these  two  calculations.  The  0.85  damping  exhibits 
slightly  less  smearing  at  the  wavefront.  In  retrospect,  very 
little  damping  would  have  been  required  in  the  linear  elastic 
analysis.  However,  hysteretic  materials  normally  exhibit  more 
severe  oscillations  due  to  changes  in  material  properties.  In 
such  calculations  heavier  damping  may  be  preferable. 


Inf luence  of  Mass  Discret i zat ion  Technique 


The  effect  of  using  the  lumped  mass  technique  in  a  two 
phase  triangular  loading  calculation  is  demonstrated  in  Figure 
3.12,  The  results  are  similar  to  those  for  the  Heaviside  loading 
which  were  discussed  in  detail  in  conjunction  with  Figure  3.4. 

Use  of  the  lumped  mass  technique  causes  heavy  smearing  in  the 
vicinity  of  the  wavefront  and  is  not  suitable  for  use  in  these 
two  phase  calculations  (refer  to  the  previous  discussion  for  more 
details ) . 
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Influence  of  Surface  Drainage  Conditions 


The  sixth  and  final  parameter  investigated  using  the 
triangular  loading  function  is  the  influence  of  the  surface 
drainage  conditions  on  the  stress  profiles.  In  the  first  portion 
of  this  phase  of  the  study,  the  triangular  loading  function  was 
applied  with  no  drainage  allowed  at  the  ground  surface.  This 
calculation  is  compared  to  th_  standard  case  in  Figure  3.13  where 
drainage  is  only  partially  inhibited  by  the  specified  pore 
pressure  loading  (i.e.  the  pore  pressure  equals  35%  of  the  total 
applied  load).  The  profiles  for  the  calculations  with  no  surface 
drainage  are  plotted  as  dashed  lines.  This  calculation  was 
performed  by  prohibiting  relative  motion  between  the  pore  water 
and  skeleton  at  the  loading  surface  while  applying  the  total  load 
to  the  ground  surface.  At  the  ground  surface,  the  pore  pressure 
at  the  impermeable  boundary  is  much  higher  than  that  on  the 
partially  drained  boundary.  Though  the  pore  pressure  builds  up 
very  rapidly  as  the  wavefront  moves  away  from  the  partially 
drained  boundary,  the  buildup  does  not  quite  reach  the  magnitude 
of  the  pore  pressure  developed  beneath  t-.e  impermeable  boundary. 
Thus,  the  pore  pressures  near  the  peak  are  slightly  higher  in  the 
case  of  the  impermeable  loading  boundary.  Because  of  the 
impermeable  boundary,  there  is  no  pore  pressure  dissipation  in 
the  near  surface  region  in  that  calculation. 
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Comparison  of  the  two  phase  calculation  on  the  impermeable 
loading  boundary  of  Figure  3.13  with  the  totally  undrained  one 
phase  calculation  of  Figure  3.7  shows  that  in  the  near  surface 
region  the  pore  pressure  response  in  both  calculations  is  nearly 
identical.  At  greater  depths,  however,  the  response  of  the  two 
phase  calculation  near  the  wavefront  matches  that  of  the  drained 
two  phase  calculations. 


Figure  3.14  shows  a  comparison  between  the  pore  pressure 
and  effective  stress  profiles  for  partially  drained  and  fully 
drained  conditions  at  the  loading  boundary.  In  the  partially 
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drained  standard  calculation,  35%  of  the  total  load  is  applied  to 
the  pore  water.  In  the  fully  drained  calculations,  the  total 
load  is  applied  on  the  soil  skeleton,  thus  maintaining  zero  pore 
pressure  at  the  loading  boundary.  In  both  calculations,  the  pore 
pressures  increase  very  rapidly  with  depth  and  are  nearly  equal 
at  the  wavefront  at  all  the  times  shown. 

It  is  concluded  from  the  results  of  the  calculations 
shown  in  Figures  3.13  and  3.14  that  wave  propagation  in  two  phase 
media  is  quite  independent  of  the  loading  and/or  drainage 
conditions  assumed  at  the  ground  surface.  This  is  a  fortunate 
effect  because  the  actual  partitioning  of  the  loading  pressure 
between  the  pore  water  and  the  soil  skeleton  on  the  surface  is 
difficult  to  assess. 


POLYNOMIAL  LOADING 


Two  parametric  calculations  were  performed  using  the 
linear  elastic  properties  of  Table  3.1  and  the  simple  fifth  order 
polynomial  pressure  loading  shown  in  Figure  3.1  The  loading  has 
a  0.1  msec  risetime  and  a  decay  described  by 


p  ( t) 


) 


(3-13) 


where  the  peak  pressure,  pQ,  is  5  ksi,  t  is  the  time,  tr  the  rise 
time,  the  decay  time,  t  ,  is  20  msec,  and  the  decay  exponent,  m, 
is  5.  This  loading  function  closely  approximates  many  explosively 
generated  airblast  loadings  with  their  rapid  early  time  pressure 
decay.  As  noted  previously,  the  impulse  under  the  fifth  order 
loading  is  i/3  of  the  impulse  under  the  triangular  loading  used 
in  the  previous  set  of  parametric  calculations. 
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One  Phase  vs.  Two  Phase  Ca leu lat ions 


Results  of  one  and  two  phase  calculations  using  the  fifth 
order  polynomial  loading  are  compared  in  Figure  3.15.  As  in  the 
previous  studies,  the  one  phase  total  stress  is  multiplied  by  0.892 
from  Equation  3.10  to  obtain  the  pore  pressure  under  the 
assumption  of  no  drainage.  The  two  phase  calculation  used  the 
standard  load  partitioning,  with  35%  of  the  total  pressure  <1750 
psi)  applied  to  the  pore  water.  Attenuation  of  the  sharp  peak  in 
the  polynomial  pressure  loading  function  is  more  rapid  than  that 
observed  for  the  triangular  loading.  At  14  msec  the  peak  pore 
pressure  from  the  one  phase  calculation  has  decayed  from  4460  psi 
at  the  surface  to  3180  psi  at  the  70  ft  depth,  a  decrease  of 
nearly  29%.  This  is  much  larger  than  the  9%  drop  observed  in  the 
corresponding  triangular  loading  calculation.  The  greater 
attenuation  is  due  to  numerical  damping  of  the  much  steeper 
pressure  decay  from  the  peak.  rJse  of  less  y  damping  would  have 
minimized  this  attenuation. 

The  two  phase  pore  pressure  profile,  shown  as  a  solid 
line  in  Figure  3,15,  is  significantly  more  smeared  and  attenuated 
than  the  one  phase  pore  pressure  profile.  This  is  similar  to  the 
effects  observed  in  both  the  Heaviside  and  triangular  loading 
cases.  The  peak  pore  pressure  in  the  two  phase  calculation  at  14 
msec  is  only  2500  psi,  44%  less  than  the  surface  peak  pore 
pressure.  The  additional  15%  drop  in  pressure  is  due  to  the 
fluid  damping  associated  with  the  relative  motion  between  the 
pore  water  and  the  soil  skeleton.  Tr.  the  triangular  loading 
calculations,  the  additional  attenuation  due  to  the  fluid  damping 
was  9%.  The  difference  indicates  that  the  fluid  damping  is 
sensitive  to  the  frequency  of  the  loading  function,  being  more 
severe  for  steeper  pressure  gradients. 
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Polynomial  vs .  Trianqu lar  Loadings 


Figure  3.16  compares  the  pore  pressure  and  effective 
stress  profiles  from  the  standard  two  phase  calculations  under 
the  polynomial  and  triangular  loadings.  The  initial  rises  in 
pore  pressure  and  effective  stress  are  identical  in  the  two 
calculations.  However,  the  more  severe  numerical  and  fluid 
damping  cause  both  the  pore  pressure  and  effective  stress  under 
the  polynomial  loading  to  attenuate  more  rapidly  than  those  under 
the  triangular  loading.  The  shape  of  the  loading  functions  and 
the  attenuation  characteristics  result  in  an  apparently  more 
rapid  propagation  of  peak  stress  under  the  polynomial  loading 
than  under  the  triangular  loading.  At  14  msec  the  peak  stress  is 
at  a  depth  of  70  ft  under  the  polynomial  loading  but  has  reached 
a  depth  of  only  63  ft  under  the  triangular  loading. 
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Table  3.1.  Material  parameters  used  in  linear  elastic 
calculations . 

MATERIAL  PROPERTIES 
(typical  values) 


Porosity 

n  ■  0.35 

Specific  gravity  of  solid  grains  (quartz) 

G  =  2.67 
s 

Permeability 

k  =  0.1  in /sec 


Unit  Weight 

Pore  water  Y  *=  62.  A  lb/ft^ 

w  3 

Bulk  mixture  y^  =  130.1  lb/ft 


Modulus 


Skeleton  bulk  modulus 

Skeleton  shear  modulus 

Skeleton  constrained  modulus 

Skeleton  Poisson's  ratio 

Water  bulk  modulus 

Solid  grain  bulk  modulus 

Mixture  undrained  bulk 
modulus  (Wood  equation) 


K 

s 

G 

s 

M 


s 


V 

s 


K 

w 


K 

m 


Decoupled  un drained  bulk 

modulus  K 

i 

Decoupled  undrained 

constrained  modulus  M 


50,111  psi 
30,067  psi 
90,200  psi 


0.25 

0.29  x  10^  psi 
5  x  10^  psi 

K  ♦  A"--—  0. 748  X  10**  psi 

W  gw 


K  +  K  -  0.798  x  10  psi 
m  s 


K  +  M  -  0.838  x  10  psi 
ms  r 
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Table  3.1.  Material  parameters  used  in  linear  elastic 
calculations  (concluded)  . 


P-wave  Speed 


Water 

c  * 
pw 

4641  ft/sec 

Skele ton 

c  = 
ps 

1583  ft/sec 

Mixture  (undrained) 

c  = 
pm 

nr 

"  =  5161  ft/sec 

Decoupled  (undrained) 

c  ,  « 
pd 

nr 

y  —  *  5463  ft/sec 
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Figure  3.3.  Pore  pressure  profiles  from  one  phase  calculations  using  consistent  and 
lumped  mass  matrices 


Figure  3.4.  Stress  profiles  from  two  phase  calculations  using  consistent  and  lumped  mass 
matrices,  Heaviside  loading 


SUI 


profiles  from  one  and  two  phase  calculations 


ms 


Figure  3.9.  Influence  of  mass  increment  factor,  r,  on  stress  profile 


r>£  numerical  y  damping  on  stress  profiles 


Sul 


using  consistent  and  lumped  mass 


Figure  3.16.  Comparison  of  stress  profiles  from  triangular  loading  and  5th  order 
polynomial  loading 


TWO  PHASE  DYNAMIC  ANALYSIS 


SATURATED  POROUS  INELASTIC  SKELETON 


INTRODUCTION 

Rischbieter  et  al.,  (1977)  and  other  investigators 
hypothesized  that  liquefaction  could  result  from  one  dimensional 
loadings  of  saturated  soils  having  hysteretic  stress  strain 
response.  An  analytic  formulation  of  this  hypothesis  was 
developed  by  Hlouin  and  Shinn  (1983)  which  illustrates  the 
liquefaction  mechanisms  using  simple  hysteretic  material  models. 
The  TP DAP  code  combines  these  mechanisms  with  the  capability  of 
handling  realistic  material  properties  to  model  wave  propagation 
in  continuous  two  phase  media.  TPDAP  permits  the  study  of 
liquefaction  using  realistic  loadings  and  in  situ  material 
profiles. 

In  this  section,  a  brief  description  of  the  liquefaction 
mechanisms  discussed  above  is  presented,  followed  by  application 
of  TPDAP  in  preliminary  studies  of  blast  induced  liquefaction  in 
two  material  types.  The  first  is  a  simple  bilinear  hysteretic 
material  model,  and  the  second  is  based  on  data  from  uniaxial 
tests  on  Enewetak  sand  presented  by  Blouin,  Martin  and  McIntosh 
(1984)  and  reported  as  a  second  volume  of  this  study. 


LIQUEFACTION  mechanisms 


The  analytic  two  phase  model  developed  by  Filouin  and 
Shinn  (1983)  treats  the  stresses  and  deformations  in  the  soil 
skeleton  as  separate  entities  during  compression  wave  loadings 
and  unloadings.  Since  granular  soil  skeletons  are  hysteretic, 
expansion  of  the  skeleton  during  unloading  will  generally  involve 
only  a  portion  of  the  volume  change  which  occurred  during  the 
compressive  wave  loadings.  During  both  loading  and  unloading, 
however,  the  pore  water  behaves  elastically.  On  unloading,  the 
pore  water  may  continue  to  expand  elastically  after  the  skeleton 
has  reached  zero  effective  stress.  At  this  point  the  expanding 
pore  water  will  begin  to  carry  the  soil  particles  into  suspension 
and  the  material  will  be  in  a  state  of  liquefaction.  Following 
liquefaction,  the  soil  will  consolidate  into  a  more  dense  state 
than  originally,  with  the  excess  pore  water  pushed  upward  into 
overlying  unsaturated  layers  or  onto  the  ground  surface.  The 
rate  of  consolidation  will  be  governed  by  the  amount  of  excess 
pore  water  (i.e.  the  extent  of  the  zone  of  liquefaction  and  the 
degree  of  liquefaction)  and  by  the  permeability  of  the  liquefied 
material  itself  and  the  permeability  of  any  overlying  non- 
liquefied  material. 

A  crude  model  depicting  both  the  soil  skeleton  and  the 

solid/water  mixture  during  uniaxial  strain  loading  and  unloading 

is  shown  in  Figure  4.1.  Loading  of  the  soil  skeleton  occurs 

along  a  bilinear  path  and  unloading  along  a  steeper  linear  path. 

The  initial  portion  of  the  loading  curve,  represented  by  the 

loading  modulus  M^,  is  assumed  to  be  elastic  and  is  limited  by  the 

maximum  strain  or  the  corresponding  effective  stress  o'  .  In 

ge  ge 

instances  where  the  elastic  limit  is  not  exceeded,  the  skeleton 
will  unload  along  the  loading  path  and  there  will  be  no  inelastic 
volume  decrease  in  the  soil  skeleton.  The  elastic  portion  of 
the  loading  curve  can  represent  either  a  cemented  soil  or  simply 
a  short  elastic  portion  of  the  loading  curve  commonly  evidenced 
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(luring  field  tests  known  as  the  seismic  or  elastic  toe.  Beyond 

the  elastic  strain  limit  e  ,  the  skeleton  is  no  longer  elastic  and 

ge  a 

loads  at  a  modulus  M 2>  and  unloads  at  a  much  stiffer  modulus  Mu* 

According  to  Terzayhi's  effective  stress  law,  the  total 

stress,  o,  applied  to  an  increment  of  saturated  soil  will  be 

carried  by  the  pore  water  pressure,  tt  ,  and  by  the  intergranular 

stress,  o'  ,  within  the  soil  skeleton.  The  pore  water  pressures 
y 

may  exceed  the  intergranular  stresses  by  an  order  of  magnitude. 

In  order  tor  liquefaction  to  occur  upon  unloading,  two  conditions 

must  be  met.  The  first  is  simply  that  the  elastic  strain  limit 

in  the  soil  skeleton  must  be  exceeded  during  loading  to  the  peak 

dynamically  applied  stress  o  .  Otherwise,  both  the  soil  skeleton 

and  pore  water  will  unload  elastically  and  no  excess  pore 

pressures  will  be  generated.  Figure  4.1  depicts  deformation  in 

the  skeleton  and  soil-water  mixture  at  depth  d  where  the  initial 

effective  overburden  stress  is  a  and  the  initial  pore  pressure 

is  iTj .  Note  that  the  in  situ  strain  in  each  case  is  taken  as 

zero  and  that  stresses  are  measured  from  the  in  situ  conditions. 

The  first  condition  for  liquefaction  is  satisfied  when  the  strain 

in  the  soil  skeleton  exceeds  the  elastic  limit  e„„.  If  no 

ge 

relative  flow  between  the  pore  water  and  skeleton  is  assumed, 

then  the  strain  in  the  solid/pore  water  mixture  equals  that  in 

the  skeleton.  Upon  unloading,  the  stiff  skeleton  unloading 

modulus  results  in  a  rapid  drop  in  effective  stress  until  it 

reaches  zero  at  the  strain  tgQ.  At  this  point,  the  skeleton 

has  lost  all  strength  and  reaches  the  liquefied  state.  Unloading 

continues  until  the  stress  in  the  soil-water  mixture  reaches  the 

effective  overburden  stress  o'  A  at  the  strain  em„.  Between 

g  u  mr 

strains  c  gQ  and  emr,  the  expanding  pore  water  tends  to  carry  the 
soil  particles  into  suspension.  Finally,  as  the  particles 
reconsolidate,  the  pore  pressure  gradually  drops  from 
its  original  value  tt.. 


back  to 


BILINEAR  MODEL 


1 


The  triangular  loading  pulse  from  Figure  3.1  was  applied 
to  a  saturated  material  having  bilinear  load-unload 
characteristics.  All  properties  of  the  skeleton  are  the  same  as 
those  used  in  the  previous  section,  given  in  Table  3.1,  except 
the  unloading  moduli  are  taken  as  three  times  the  corresponding 
loading  moduli.  Thus,  hysteresis  is  introduced  into  the  material 
skeleton  with  a  strain  recovery  ratio  of  1/3. 

The  results  of  this  calculation  are  plotted  as  a  solid 
line  in  Figure  4.2.  The  effective  stress  and  pore  pressure 
profiles  at  four  times,  from  2  to  14  msec,  are  compared  to  those 
from  the  corresponding  calculations  using  the  linear  elastic 
model.  At  2  msec  there  is  no  apparent  difference  between  the 
calculations .  At  later  times,  however,  the  different  unloading 
slopes  affect  the  response  behind  the  stress  peaks.  This  effect 
becomes  more  pronounced  with  increasing  depth  and  time.  As 
explained  in  the  previous  subsection,  during  unloading  the 
bilinear  skeleton  will  tend  to  recover  only  1/3  of  the  loading 
strain,  and  the  skeleton  stress  will  drop  proportionately  faster. 
For  a  totally  undrained  loading  using  this  set  of  bilinear 
properties,  the  effective  stress  would  be  expected  to  drop  to 
zero  about  1/3  of  the  way  into  the  unloading.  At  14  msec  the 
effective  stress  in  Figure  4.2  drops  to  zero  at  a  depth  of  about 
22  ft,  at  approximately  1/3  of  the  way  into  the  unloading.  At 
this  point,  the  soil  particles  would  tend  to  begin  separating  and 
a  state  of  liquefaction  would  exist.  In  this  bilinear 
calculation,  however,  tension  is  allowed  to  develop  in  the  soil 
skeleton,  so  the  post-liquefaction  process  is  not  accurately 
depicted. 
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Because  of  the  hysteresis  in  the  bilinear  skeleton,  the 
effective  stress  drops  faster  during  unloading  than  the  effective 
stress  in  the  elastic  case.  The  pore  pressure  behind  the 
wavefront  in  the  bilinear  material  drops  less  rapidly  than  the 
pore  pressures  in  the  linear  material,  compensating  for  the  more 
rapid  effective  stress  drop  in  the  bilinear  soil.  The  total 
stresses  are  about  equal  in  both  calculations. 


.  -< 


The  inelastic  behavior  of  the  skeleton  also  causes 
additional  energy  attenuation  due  to  the  hysteretic  material 
damping.  Careful  examination  of  the  14  msec  profile  in  Figure 
4.2  shows  that  the  total  stress  behind  the  wavefront  is  lower  in 
the  bilinear  material  than  that  in  the  elastic  case.  This 
difference  evidently  results  from  the  hysteretic  damping. 


Figure  4.3  is  an  expanded  view  of  the  bilinear  effective 
stress  and  pore  pressure  profiles  at  15,  16  and  17  msec.  A  zone 
of  effective  tensile  stress  is  developing  at  these  later  times. 
This  zone  extends  from  a  depth  of  about  14  ft  to  a  depth  of  40  ft 
at  a  time  of  17  msec.  This  developing  tensile  zone  is  indicative 
of  liquefaction,  but  because  tensile  stresses  are  allowed  to 
develop  in  the  skeleton  and  because  there  are  no  gravitational 
stresses  in  this  calculation,  it  does  not  realistically  model  the 
liquefaction  process. 


E NEW ETA K  SAND 

The  final  set  of  parametric  calculations  undertaken  in 
this  study  examines  the  response  of  saturated  sand  from  Enewetak 
Atoll.  The  skeleton  properties  are  taken  from  laboratory  data 
reported  in  the  second  volume  of  this  study  by  Blouin,  Martin  and 
McIntosh  (1984).  The  sand  is  a  uniform  carbonate  beach  sand  with 
a  mean  grain  diameter  of  about  0.5  mm  and  a  porosity  of  45%. 


Properties  selected  for  this  sand  are  shown  in  Figure  4.4.  The 
loading  modulus  is  a  stepwise  approximation  to  the  actual 
constrained  compression  data  from  test  12.  The  unloading  modulus 
is  expressed  as  a  function  of  peak  stress  according  to 
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u 


3470  a 


(0.651) 

am 


(4-1) 


where  o'  is  the  maximum  previous  axial  stress  expressed  in  psi. 

am  ^ 

As  indicated  by  the  data,  the  unload  modulus  becomes  continually 
stiffer  with  increasing  peak  stress. 

Three  uniaxial  strain  calculations  were  performed  on  the 
saturated  Enewetak  sand.  All  used  a  triangular  loading  with 
duration  and  risetime  as  specified  in  Figure  3.1.  The  peak 
stress  in  the  first  two  calculations  was  50  ksi  and  the  peak 
stress  in  the  last  calculation  was  5  ksi.  Unlike  the  bilinear 
calculations,  tensile  stress  was  not  allowed  to  develop  in  the 
Enewetak  sand.  Whenever  the  effective  stress  dropped  to  zero, 
the  skeleton  had  no  strength  and  the  modulus  dropped  to  zero. 


50  ks i  Peak  Stress 

The  two  50  ksi  calculations  examined  the  influence  of  the 
surface  drainage  conditions.  Neither  calculation  included  in 
situ  stresses  due  to  gravity.  In  addition  to  the  usual  effective 
stress  and  pore  pressure  profiles,  motion  and  stress  time 
histories  at  selected  depths  are  also  presented.  Figure  4.5 
shows  the  stress  profiles  from  the  two  calculations  at  2,  6,  10 
and  14  msec.  Because  the  skeleton  modulus  of  the  Enewetak  sand 
was  considerably  lower  than  that  of  the  previous  skeletons 
analyzed  in  this  study,  the  effective  stress  profiles  are 
expanded  by  a  factor  of  ten  in  Figure  4.5  to  afford  better 
definition.  The  most  significant  feature  of  both  sets  of 


effective  stress  profiles  is  the  development  of  a  large  zone  of 
liquefaction  behind  the  wavefront.  By  14  msec  this  zone  extends 
from  the  surface  to  a  depth  of  42  ft. 

There  is  little  difference  between  the  profiles  for  the 
two  different  surface  drainage  conditions.  The  partial  surface 
drainage  condition  was  attained  by  partitioning  the  airblast 
loading  between  the  pore  water  and  soil  skeleton  in  proportion  to 
the  porosity.  Thus,  45%  of  the  total  load  was  applied  to  the 
pore  water  and  the  remaining  portion  to  the  soil  skeleton.  In 
the  calculation  with  no  surface  drainage,  relative 
motion  between  the  soil  skeleton  and  pore  water  was  prohibited  at 
the  loading  boundary.  As  was  also  observed  in  the  linear 
elastic  analysis,  the  surface  stress  distribution  is  rapidly 
repartitioned  so  that  the  undrained  surface  loading  profile  is 
approximated  within  a  very  short  distance  beneath  the  surface. 

In  other  words,  the  surface  load  partitioning/drainage  conditions 
affect  only  the  very  near  surface  response.  The  overall  response 
is  insensitive  to  the  assumed  surface  drainage  conditions. 

Time  histories  of  the  response  of  the  soil  skeleton  at  a 
depth  of  10  ft  are  shown  for  both  the  partial  surface  drainagi 
and  no  surface  drainage  calculations  in  Figure  4.6.  The  top  time 
history  compares  the  effective  vertical  and  horizontal  stresses 
from  the  two  calculations,  and  the  remaining  time  histories 
compare  the  displacement,  velocities  and  accelerations.  As 
mentioned  in  the  above  discussion  of  the  stress  profiles,  the 
influence  of  surface  drainage  is  small.  According  to  the 
effective  stress  time  histories,  liquefaction  occurs  at  the  10  ft 
depth  shortly  after  6  msec.  The  peak  vertical  effective  stress 
in  the  skeleton  is  approximately  1400  psi  and  the  peak  horizontal 
effective  stress,  controlled  by  the  K  value  of  0.5,  is  about  700 
psi.  Thus,  the  peak  skeleton  stress  at  this  depth  is  only  2.8% 
of  the  50  ksi  peak  loading  stress.  Since  so  little  stress  is 
carried  in  the  skeleton,  the  influence  of  the  liquefaction  on  the 
overall  skeleton  motions  during  the  constrained  dynamic  loading 
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is  negligible.  However,  it  should  be  noted  that  the  material  in 
the  liquefied  condition  no  longer  has  any  shear  strength.  This 
can  result  in  very  large  mass  movements  at  late  times  under 
gravitational  or  other  loadings  whenever  nonuniform  stress  fields 

exist. 


The  corresponding  pore  pressure  time  histories  and 
apparent  relative  fluid  motions  between  the  pore  water  and  soil 
skeleton  at  the  10  ft  depth  are  plotted  in  Figure  4.7.  As 
described  in  the  theoretical  discussion  of  Section  2,  the 
apparent  relative  fluid  motion  equals  r.  times  the  actual  relative 
motion  where  n  is  the  porosity.  In  the  Enewetak  sand,  the  actual 
relative  motions  between  the  particles  and  pore  water  is  2.2 
times  the  apparent  values.  The  pore  pressure  time  histories 
resemble  the  triangular  loading  pulse.  Peak  pressures  are  about 
43  ksi,  over  30  times  higher  than  the  peak  vertical  stress  in  the 
soil  skeleton  at  this  depth.  The  displacement  time  history 
indicates  that  the  pore  water  moves  downward  relative  to  the  soil 
skeleton  during  the  dynamic  loading.  The  peak  magnitude  of 
apparent  relative  displacement  is  approximately  0.9  in  and  occurs 
just  after  6  msec,  which  also  is  the  time  at  which  liquefaction 
occurs . 


The  velocity  time  histories  indicate  a  sharp  downward 
velocity  relative  to  the  soil  skeleton  at  the  wavefront.  The 
peak  apparent  relative  downwarc  velocity  is  about  47  ft/sec  and 
the  actual  peak  relative  velocity  is  about  103  ft/sec.  This 
occurs  at  2.4  msec.  It  compares  to  a  peak  downward  velocity  in 
the  soil  skeleton  of  approximately  325  f.  t/sec,  though  this  peak 
is  reached  somewhat  later  at  about  4  msec.  The  apparent  relative 
velocity  decays  rapidly  from  its  peak  value  and  levels  off  at 
about  4  ft/sec  by  5  msec.  Just  aftec  6  msec,  at  the  time  of 
liquefaction,  the  velocity  suddenly  reverses  direction  and 
becomes  stable  in  the  upward  direction  with  a  magnitude  of  5 
ft/sec.  This  reversal  is  evidently  associated  with  the 
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liquefaction.  The  reversal  in  velocity  controls  the  time  of  peak 
in  the  apparent  relative  displacement  time  history. 

The  apparent  relative  acceleration  is  initially  downward 
with  a  peak  value  of  about  2700  gs.  This  compares  to  a  peak 
skeleton  acceleration  of  14000  gs.  The  apparent  relative 
acceleration  quickly  reverses  at  2.4  msec  and  reaches  a  peak 
deceleration  of  about  900  gs  at  2.8  msec.  It  then  tapers  off  to 
near  zero  slightly  after  5  msec  and  registers  a  small  negative 
value  at  the  time  of  liquefaction. 

Figures  4.8  and  4.9  present  stress  and  motion  time 
histories  for  the  50  ksi  triangular  loading  at  a  depth  of  20  ft 
for  the  undrained  surface  loading  condition.  These  waveforms  are 
similar  to  those  at  the  10  ft  depth  of  Figures  4.6  and  4.7.  As 
shown  in  the  effective  stress  time  histories  of  Figure  4.8, 
liquefaction  occurs  at  8.4  msec  at  the  20  ft  depth.  The  peak 
apparent  relative  velocity  between  the  pore  water  and  the 
skeleton  is  41  ft/sec  at  this  depth  and  the  peak  apparent 
relative  acceleration  is  about  1700  gs.  As  was  the  case  at  the 
10  ft  depth,  the  relative  velocity  reverses  at  the  time  of 
liquefaction  indicating  that  the  pore  water  begins  to  move  upward 
relative  to  the  soil  skeleton. 

_5  ksi  Peak  St ress 


The  effective  stress  and  pore  pressure  profiles  as 
functions  of  time  and  depth  for  the  final  Enewetak  sand 
calculation  are  presented  in  Figure  4.10.  The  applied  loading 
was  the  5  ksi  triangular  pulse  from  Figure  3.1  with  no  surface 
drainage  permitted  at  the  loading  boundary.  In  this  calculation 
gravitational  stresses  were  included  on  both  the  soil  skeleton 
and  the  pore  water.  The  initial  in  situ  effective  stress  is 
apparent  in  all  the  profiles,  where  the  effective  stress 
increases  linearly  with  depth  and  at  100  ft  equals  43  ps i . 
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Coincidentally,  the  pore  pressure  under  the  gravitational  loading 
also  is  43  psi  at  that  depth.  Because  the  effective  stress 
profiles  are  magnified  by  a  factor  of  10  compared  to  the  pore 
water  pressures,  the  43  psi  pore  pressure  is  barely  noticeable  in 
the  scale  of  Figure  4.10. 

Several  differences  are  apparent  between  the  5  ksi 
loading  profiles  of  Figure  4.10  and  the  50  ksi  loading  profiles 
of  Figure  4.5.  At  10  msec  the  peak  effective  stress  in  the  5  ksi 
case  is  only  about  1.6%  of  the  peak  pore  pressure.  In  the  50  ksi 
loading  at  this  time,  the  peak  effective  stress  is  about  3.6%  of 
the  peak  pore  pressure.  The  significantly  greater  portion  of  the 
total  stress  carried  in  the  pore  skeleton  under  the  50  ksi 
loading  is  due  to  the  increase  in  skeleton  modulus  with 
increasing  strain  indicated  in  the  skeleton  properties  of  Figure 
4.4.  Under  the  50  ksi  loading  the  skeleton  strains  are 
proportionally  higher  than  those  under  the  5  ksi  loading;  thus, 
the  increasingly  stiff  skeleton  modulus  at  these  higher  strains 
results  in  a  relatively  greater  portion  of  the  total  stress 
distributed  to  the  skeleton. 

Another  difference  observed  in  Figure  4,10  is  the 
increase  in  dynamic  peak  effective  stress  with  increasing  time 
and  depth.  This  is  in  contrast  to  the  decreasing  peak  effective 
stress  in  the  50  ksi  loading.  Here  the  difference  is  due  to  the 
application  of  the  gravitational  effective  stress  in  the  5  ksi 
case.  With  increasing  depth,  the  initial  in  situ  effective 
stress  increases  and  the  skeleton,  in  essence,  becomes  stiffer 
with  increasing  depth.  Thus,  for  a  given  peak  total  stress 
level,  the  portion  of  the  total  stress  carried  in  the  skeleton 
would  be  less  at  shallow  depths  than  it  would  be  at  deeper  depths. 
However,  because  peak  stress  attenuates  with  increasing  depth  due 
to  hysteretic  and  fluid  damping,  the  net  change  in  peak  effective 
stress  depends  on  the  summation  of  these  two  effects.  In  this 
case,  the  net  result  is  an  increase  in  effective  stress  with 
increasing  depth. 
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In  the  5  ksi  loading  calculations  with  gravitational 
stresses,  the  soil  also  liquefied  during  unloading  from  the  peak 
dynamic  stress.  As  shown  in  Figure  4.10,  at  14  msec  the  soil  is 
liquefied  from  the  surface  to  a  depth  of  shout  36  ft.  This  is 
somewhat  less  than  the  42  ft  depth  of  liquefaction  shown  in 
Figure  4.5  for  the  50  ksi  load  at  14  msec. 

I 

Effect-  stress  and  skeleton  motion  time  histories  at 
the  10  ft  and  20  ft  depths  for  the  5  ksi  loading  are  shown  in 
Figure  4,11  and  4.13  respectively.  The  pore  pressure  and 
apparent  relative  motion  time  histories  for  the  10  ft  and  20  ft 
depths  are  shown  in  Figures  4.12  and  4.14.  The  waveforms  in 
these  examples  are  similar  to  those  presented  ir.  the 
corresponding  Figures  4.6  througn  4.9  for  me  50  ksi  loading. 

The  magnitudes  are  correspondingly  lower  for  the  5  ksi  loading. 
The  only  significant  difference  appears  to  oe  that  the  apparent 
relative  velocity  of  the  pore  water  becomes  negative  relative  to 
the  skeleton  prior  to  liquefaction  in  the  5  ksi  case.  Thus, 
velocity  reverses  gradually  during  the  initial  unloading  rather 
than  abruptly  at  the  time  of  liquefaction. 
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Figure  4.4.  Loading  and  unloading  skeleton  constrained  moduli  of  Enewetak  sand 
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Figure  4.10.  Stress  profiles  under  5  ksi  triangular  loading  of  Enewetak  sand 
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Figure  4.12. 


Pore  pressure  and  apparent  relative  fluid  motion 
time  histories  at  10  ft  depth,  5  ksi  triangular 
loading 
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SECTION  5 


SUMMARY  AND  CONCLUSIONS 


The  dynamic  response  of  saturated  porous  media  is 
examined  under  uniaxial  loadings.  The  material  models  used  in 
this  study  become  increasingly  more  sophisticated  culminating 
with  an  actual  saturated  sand  from  Enewetak  Atoll.  In  Section  2, 
the  theoretical  background  and  numerical  code,  TPDAP ,  used  in 
this  study  are  described.  The  governing  differential  equations 
of  motion  for  the  bulk  soil-water  mixture  and  the  pore  fluid  are 
derived.  Hiot's  fundamental  analytic  work  describing  the 
behavior  of  saturated  porous  media  is  reviewed  and  summarized. 
His  differential  equation  for  relative  fluid  motion  is 
generalized  by  inclusion  of  the  mass  increment  factor  r  of 
Equation  2-60  which  accounts  for  additional  friction  force 
resulting  from  relative  acceleration  between  the  pore  water  and 
soil  skeleton.  ttiot  derived  values  of  r  for  two  simple  assumed 
pore  geometries  but  actual  values  of  r  for  real  soils  have  not 
been  determined.  Values  of  r  for  soils  should  be  obtained 
through  a  combination  of  further  analytical  studies  and 
laboratory  investigations.  The  latter  portion  of  Section  2 
presents  the  finite  element  formulations  used  in  the  TPDAP  two 
phase  dynamic  analysis  program. 

Section  d  presents  the  results  of  dynamic  uniaxial 
loadings  of  saturated  porous  media  having  a  linear  elastic 
skeleton.  An  analysis  of  these  results  produced  the  following 
conclu  s  ions : 
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Use  of  lumped  mass  discretization  techniques  are  not 
suitable  tor  two  phase  dynamic  analysis.  Instead, 
consistent  mass  techniques  are  recommended; 

Surface  drainage  conditions  and/or  partitioning 
of  the  applied  load  between  the  pore  fluid  and 
soil  skeleton  only  influence  the  dynamic  response 
in  the  immediate  vicinity  of  the  loaded  boundary. 
Below  this  zone,  the  dynamic  response  is 
essentially  independent  of  the  assumed  loading 
distribution  and  boundary  drainage  conditions.  To 
avoid  undesirable  numerical  oscillations  caused 
by  large  stress  gradients,  use  of  an  impermeable 
loading  boundary  is  generally  recommended; 


Relative  motion  between  the  pore  fluid  and  soil 
skeleton  near  the  wavefront  creates  fluid 
frictional  dissipation  which  tends  to  smear  the 
wavefront  and  attenuate  the  peak  stress.  These 
effects  are  a  function  of  permeability  and  pore 
geometry  as  represented  by  the  mass  increment 
factor,  r.  As  permeability  decreases,  the  dynamic 
response  approaches  that  obtained  using  an  equiv¬ 
alent  one  phase  analysis.  As  r  increases  the  rel¬ 
ative  motions  decrease  and  the  dynamic  response 
approaches  that  of  an  equivalent  one  phase  loading. 
Values  of  r  of  1/3  and  1/5,  derived  for  the 
simplified  circular  and  flat  pore  geometries, 
had  a  negligible  influence  on  the  overall  dynamic 
response  in  these  calculations.  Larger  values 
of  r  which  may  be  representative  of  real  soils 
could  have  a  significant  influence  on  overall 
response ; 


•  fluid  damping  results  in  more  rapid  stress 
attenuation  under  loadings  with  more  rapid 
pressure  decays; 

•  Artificial  numerical  y  damping  eliminated  high 
frequency  oscillations,  but  smeared  the  wavefront 
at  the  toe.  Minimal  artificial  damping  would  have 
been  required  in  these  elastic  ca  l'-n  lat  ions . 

In  Section  4,  dynamic  uniaxial  loadings  were  applied  to 
saturated  porous  media  having  a  bilinear  hysteretic  skeleton  and 
the  skeleton  properties  of  Enewetak  sand.  Principal  con¬ 
clusions  and  observations  include: 

•  Liquefaction  occurred  in  the  calculations  using  both 
the  bilinear  and  Enewetak  hysteretic  soil  skeletons 
during  the  unloading  portion  of  the  dynamic  response. 
Liquefaction  resulted  from  non-recoverable  strains 
induced  in  the  skeleton  during  loading. 

•  The  TP DAP  two  phase  code  provides  a  powerful  tool  for 
studying  wave  propagation  in  two  phase  materials. 

It  permits  detailed  examination  of  dynamic  material 
response  including  pore  pressure,  effective  stresses, 
skeleton  motions  and  relative  fluid  motions  at  specified 
points  in  time  and  space.  based  on  the  relative  fluid 
motions,  pore  water  migration  and  flow  rate  can  be 
computed . 

•  Relationships  between  pore  pressure,  effective  stress, 
skeleton  motion  and  relative  motion  between  the  pore 
fluid  and  skeleton  were  examined  in  the  calculations  of 
the  Enewetak  sand.  During  loading,  motion  of  the  pore 
fluid  relative  to  the  skeleton  is  in  the  direction  of 
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the  wave  propagation.  During  the  unloading,  the  relative 
velocity  rapidly  decreases  and  has  reversed  by  the  time 
of  liquefaction  so  that  relative  motion  is  opposite  to 
the  direction  of  propagation. 

In  situ  gravitational  stresses  can  result  in  an 
increase  in  dynamic  effective  stress  during  loading 
with  increasing  depth  and  time.  The  occurrence 
of  such  an  increase  is  dependent  on  the  stress 
strain  properties  of  the  soil  skeleton  and  the 
attenuation  character istics  of  the  soil. 
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